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A COMPARATIVE STUDY OF THE HARRIS-PRIESTER , 
JACCHIA-ROBERTS, AND MS IS ATMOSPHERIC DENSITY 
MODELS IN THE CONTEXT OF SATELLITE ORBIT 
DETERMINATION* 

R. E. Shanklin, Jr., T. Lee, M. K. Mallick, R. A. Kuseski, 

and J. 0. Cappellari, Jr. 

Computer Sciences Corporation 

ABSTRACT 

Extensive comparisons of the Harris-Priester, Jacchia- 
Roberts, and MSIS (Mass Spectrometer/Incoherent Scatter) 
atmospheric density models as used in satellite orbit deter 
mination are summarized. The quantities compared include 
Bayesian weighted least squares differential correction sta 
tistics and orbit solution consistency and accuracy. 


*This work was supported by the Operations Analysis Section 
Operational Orbit Support Branch, Goddard Space Flight 
Center, National Aeronautics and Space Administration, 
under Contract NAS 5-24300. 
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SECTION 1 


INTRODUCTION 


Atmospheric drag is a significant perturbation of Earth sat- 
ellite orbits with perigee heights of less than 1000 kilom- 
eters. The acceleration of a spherical satellite due to 
atmospheric drag is given by the equation 


drag 


1 S 

2 m 


A p |V| V 


where p = atmospheric density at the position of the satel- 
lite 

V = satellite velocity relative to the atmosphere 

A = satellite reference cross-sectional area 

Cp = satellite drag coefficient 

m = satellite mass 


Therefore, calculation of the drag acceleration requires 
knowledge of the atmospheric density as a function of posi- 
tion and time . 


This paper presents the results of a comparative study of 
three different global atmospheric density models in the 
context of orbit determination. The three models compared 
are the Harr is-Pr iester (H-P) model , the Jacchia-Rober ts 
(J-R) model, and the Mass Spectrometer/Incoherent Scatter 
(MS IS) model. 

The Harr is-Pr iester model is based on theoretical tempera- 
ture profile solutions of the heat conduction equation under 
hydrostatic equilibrium conditions. The model assumes two 
heat sources; solar extreme ultraviolet (EUV) heating and 
an artificial heat source that produces the diurnal varia- 
tion deduced from satellite drag calculations. In the mod- 
ified Harr is-Pr iester model used for this study, the EUV 


1-2 



heating level is selected by choosing among 10 different 
altitude-density profile tables representing 10 different 
levels of solar , flux, and the diurnal variation is modeled 
by a correction calculated using a power of a cosine 
(References 1 and 2) . 

The Jacchia-Rober ts model is based on empirical temperature 
profiles scaled by an upper boundary exospheric temperature 
(Tj^) . Analytic density calculation is accomplished through 
integration of thermodynamic equations. The modeling in- 
cludes corrections for EUV heating, solar particle flux 
(so-called geomagnetic) heating, semiannual variations, sea- 
sonal variations, and the diurnal variation (References 2 
and 3) . 

The MSIS model is based on fitting spherical surface har- 
monic expansions to match the angular dependence exhibited 
by mass spectrometer and incoherent scatter measurements. 

The MSIS formulation includes sections that model EUV heat- 
ing, solar particle flux heating, annual variations, semian- 
nual variations, diurnal variations, semidiurnal variations, 
terdiurnal variations, and departures from diffusive equi- 
librium. MSIS modeling has been implemented in a special 
GTDS load module. Dr. Hedin and his associates at the 
Goddard Space Flight Center, who developed the model (Ref- 
erence 4) , contributed advice and some of their program sub- 
routines during the GTDS implementation. 

Table 1 shows sample density profiles for the three atmos- 
pheric models with two different solar EUV levels and one 
geomagnetic activity level. Figure 1 shows the Jacchia- 
Roberts and MSIS densities, relative to the Harr is-Pr iester 
density, as a function of altitude. The figure shows max- 
imum ratios as high as 2.0 but, as is apparent from the 
table, the three profiles are quite similar in overall shape. 
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TABLE 1. ATMOSPHERIC DENSITIES COMPUTED USING HARRIS-PRIESTER, 
JACCHIA-ROBERTS, AND MSIS MODELS 



DENSITY (kg/km^) 

ALTITUDE 

(km) 

HARRIS-PRIESTER 

JACCHIA-ROBERTS 

MSIS 

F,0.7 = 125.0 

Fio7 = 150.0 

F,o.7 = 116.2 

Fio.7 = 135.1 

Fio7 = 140.0 
F,o.7 = 165.3 

^0.7 = ’ '6-2 
Fio.7= 135.1 

Fio7 = 140.0 

Fio.7 = 165.3 

150 

.205 E + 1 

.206 E + 1 

.193 E + 1 

.210 E + 1 

.203 E + 1 

.204 E + 1 

200 

.224 E 0 

.255 E 0 

.228 E 0 

.270 E 0 

.274 E 0 

.313 E 0 

250 

.459 E - 1 

.583 E - 1 

.559 E - 1 

.721 E - 1 

.636 E - 1 

.802 E - 1 

300 

.129 E - 1 

.178 E-1 

.177 E - 1 

.249 E - 1 

.187 E-1 

.255 E - 1 

350 

.425 E - 2 

.631 E - 2 

.637 E - 2 

.977 E - 2 

.633 E - 2 

.926 E - 2 

400 

.155 E -2 

.247 E - 2 

.246 E - 2 

.413 E - 2 

.236 E - 2 

.368 E - 2 

450 

.521 E - 3 

.879 E - 3 

.835 E-3 

.157 E -2 

.780 E-3 

.131 E - 2 

500 

.218 E -3 

.392 E - 3 

.353 E-3 

.724 E-3 

.324 E-3 

.582 E-3 

550 

.963 E - 4 

.182 E-3 

.1 55 E - 3 

.344 E -3 

.139 E-3 

.266 E-3 

600 

.451 E - 4 

.851 E - 4 

.706 E - 4 

.169 E-3 

.619 E - 4 

.125 E-3 

650 

.227 E - 4 

.451 E -4 

.339 E - 4 

.851 E - 4 

.285 E - 4 

.600 E - 4 

700 

.112E-4 

.21 7 E - 4 

.154 E - 4 

.394 E - 4 

.120 E - 4 

.259 E -4 

750 

.691 E - 5 

.127 E -4 

.878 E-5 

.219 E -4 

.623 E-5 

.134 E -4 

800 

.464 E - 5 

.804 E - 5 

.548 E-5 

.128 E - 4 

.352 E-5 

.728 E-5 

850 

.316 E - 5 

.462 E - 5 

.348 E - 5 

.737 E-5 

.200 E-5 

.378 E-5 

900 

.245 E - 5 

.301 E - 5 

.258 E-5 

.500 E-5 

.137 E-5 

.236 E-5 

950 

.198E-5 

.201 E - 5 

.201 E-5 

.361 E-5 

.102 E-5 

.158 E-5 

1000 

.163 E -5 

.141 E-5 

.155 E-5 

.262 E-5 

.761 E - 6 

.107 E-5 


NOTES: 1. Kp = 3.3 FOR JACCHIA-ROBERTS DENSITY AND Ap = 33 FOR MSIS DENSITY ARE USED. 

2. THESE PROFILES ARE FOR AUGUST 30, 1978, AT A LATITUDE OF 46° N, AN EAST LONGITUDE OF 205°, 
AND A LOCAL SOLAR TIME OF 1:40 P.M. 
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FIGURE 1. ATMOSPHERIC MODEL DENSITY RATIOS 
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SECTION 2 - COMPARATIVE STUDY STRUCTURE 

All the results presented in Section 3 of this paper are 
based on Goddard Trajectory Determination System (GTDS) 
Bayesian weighted least squares differential correction so- 
lutions. Nine different series of six GTDS Differential 
Correction (DC) Program runs were made for each of the three 
atmospheric models. Three different satellites, with per- 
igee heights between 310 and 560 kilometers, were studied; 
other orbital parameters for these satellites are given in 
Table 2. The nine series of orbit determination arcs are 
listed in Table 3. 

Each series contains six 30-hour-arc solutions. The solu- 
tions are used to generate 30-hour ephemerides that overlap 
adjacent ephemerides by 6 hours. The ephemerides are then 
compared in order to determine the maximum position differ- 
ences (in the orbital reference frame) during the overlap 
periods. The 162 DC Program solutions produce 135 maximum 
overlap position differences. These differences are used to 
evaluate the consistency and accuracy obtained when each of 
the three atmospheric models is used. 

Each differential correction solution is made up of seven 
numbers: three position coordinates, three velocity coor- 

dinates, and the drag variation parameter ( p^_) , which is 
a scaling factor in the drag acceleration equation, i.e., 

I ^ »i) ^ 

This scaling factor is applied during generation of the 
ephemeris that uses the differential correction solution. 


^drag 
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TABLE 2. SATELLITE ORBtTAL ELEMENTS 


SATELLITE 

DATE 

PERIGEE 

HEIGHT 

(kilometers) 

APOGEE 

HEIGHT 

(kilometers) 

INCLINATION 

(degrees) 

AE-3 

AulsUST 1, 1978 

331 

341 

68 

MAG SAT 

OCTOBER 31, 1979 

352 

561 

97 


MARCH 1, 1980 

323 

471 

97 

SAGE 

FEBRUARY 19, 1979 

560 

655 

55 


TABLE 3. COMPARATIVE STUDY SERIES 


SERIES 

NUMBER 

SATELLITE 

TIMESPAN 

1 


AUGUST 1-6, 1978 

2 

AE-3 

AUGUST 14-19, 1978 

3 


SEPTEMBER 2-8, 1978 

4 


OCTOBER 31 -NOVEMBER 5, 1979 

5 


DECEMBER 1 -6, 1979 

6 

MAGSAT 

JANUARY 1-6, 1980 

7 


FEBRUARY 1-6, 1980 

8 


MARCH 1-6, 1980 

9 

SAGE 

FEBRUARY 19-25, 1979 
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Spacecraft attitude is not considered, since a spherical 
model is employed. Furthermore, no aerodynamic forces 
(e.g. , lift) other than drag are modeled. The spherical 
approximation is crude for all three satellites, and it is 
possible that other aerodynamic forces are nonnegligible . 
However, it is reasonable to expect that both assumptions 
have a negligible effect on the results of this study, 
because the results are obtained by applying each of the 
three atmospheric models to the same arcs with the same ob- 
servation sets. Simply stated, unmodeled aerodynamic forces 
should perturb the solutions for all three atmospheric 
models in a similar manner . 
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SECTION 3 - COMPARATIVE STUDY RESULTS 

This section summarizes the results of this comparative 
study of atmospheric density models in the context of 
short-arc (30-hour) orbit determination. A detailed, run- 
by-run presentation of these results is available in Ref- 
erence 5. Two cautionary remarks are appropriate. 

First, these results should not be interpreted as a compar- 
ison of atmospheric models; conclusions about the relative 
merits of the models must be limited to this highly spe- 
cialized context--short-arc orbit determination in which an 
average drag scaling factor is solved for. 

Second, any series of orbit determination and ephemeris com- 
parison runs may contain a few sporadic large overlap dif- 
ferences and a few differential corrections with large RMS 
residuals. Some of the runs included in this study show 
such large differences and/or high RMSs. 

The average weighted RMSs and the average maximum position 
differences for the three AE-3 series are given in Table 4. 
The averages over all three series are also given, along 
with the ranges of the EUV heating index , (Fj_q and the 
solar particle flux index (Kp) . The averages ‘show that 
the Jacchia-Rober ts overlap differences are about 11.5 per- 
cent (24 meters) smaller than the Harr is-Pr iester averages 
and that the MSIS averages are about 19 percent (38 meters) 
larger than the Harr is-Pr iester averages. The 62-meter dif- 
ference between the Jacchia-Rober ts and MSIS averages cannot 
be considered either large or significant. 

The same information is given for Magsat in Table 5. This 
study includes five series of arcs. The Magsat results show 
that both the Jacchia-Rober ts and MSIS average differences 
are about 9 percent larger than the Harr is-Pr iester average 
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TABLE 4. COMPARATIVE ATMOSPHERIC DENSITY MODEL STUDY RESULTS 
FOR AE-3 (AUGUST AND SEPTEMBER 1978) 


SERIES 

RANGE OF F^q^ 
watt/(m^ -Hz)j 

RANGE OF 

HARRIS-PRIESTER MODEL 

JACCHIA-ROBERTS MODEL 

MSIS 

VIODEL 

AVERAGE 

WEIGHTED 

RMS 

MAXIMUM 

POSITION 

DIFFERENCE 

(meters) 

AVERAGE 

WEIGHTED 

RMS 

MAXIMUM 

POSITION 

DIFFERENCE 

(meters) 

AVERAGE 

WEIGHTED 

RMS 

MAXIMUM 

POSITION 

DIFFERENCE 

(meters) 

AUGUST 1-6 

106.0-117.6 

0-6 

■■ 

191 

5.2 

175 

8.4 

265 

AUGUST 14-19 

115.6-134.9 

0-6 

mm 

225 

7.8 

217 

8.5 

324 

SEPTEMBER 2-8 

159.8-181.1 

0-6 

M 

209 

8.4 

163 

7.2 

164 

averages 

- 

- 

6.5 

208 


184 

8.0 

251 
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TABLE 5. COMPARATIVE ATMOSPHERIC DENSITY MODEL SiUDY RESULTS FOR MAGSAT 
(NOVEMBER AND DECEMBER 1979; JANUARY, FEBRUARY, AND MARCH 1980) 


FEHIOO 


RANGE 
OP PlO.7 
VARIATION 
0 0 ■ 22 walls m' ^ 1 ) 


RANGE 
OF Kp 
VARIATION 


HARRIS PRIESTER RESULTS 


WEIGHTED 

RMS 


MAXIMUM 

POSITION 

DIFFERENCE 

(meters) 


JACCHIA- ROBERTS RESULTS 


WEIGHTED 

RMS 


MAXIMUM 

POSITION 

DIFFERENCE 

(meters) 


MSIS RESULTS 


WEIGHTED 

RMS 


MAXIMUM 

POSITION 

DIFFERENCE 

(meters) 


OCT. 31 NOV. 5. 1979 
DEC. 1 -6, 1979 
JAN. 1 6. 1980 
FEB. 1- 6, 1980 
mar. 1 6. 1980 


AVERAGES 


207.5- 214.9 

152.2 223.4 
188.9 212.4 

212.6- 231.7 

170.2 176.7 


0 4 

0 -4 
1-5 
0-4 
0 3 


8.3 
12.4 

9.4 
12.7 

9.8 


204 

213 

326 

161 


7.8 

11.5 
9.5 

12.5 

13.4 


175 

166 

298 

396 


10.6 


222 


242 


8.0 

12.8 

11.3 

13.8 

10.0 


11.2 


190 

255 

288 

313 
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differences. As in the case of AE-3, the Magsat results 
demonstrate that the three atmospheric density models are 
comparable in the context of this study. 

The average RMSs and overlap position differences for the 
series of SAGE arcs are given in Table 6. Both the RMSs and 
the overlap differences agree to within 3 percent; all three 
atmospheric models produce essentially equivalent errors. 

TABLE 6. COMPARATIVE ATMOSPHERIC DENSITY MODEL STUDY RESULTS 
FOR SAGE (FEBRUARY 19-25, 1979) 


ATMOSPHERIC 

DENSITY 

MODEL 

USED 

AVERAGE 

WEIGHTED 

RMS 

AVERAGE 

MAXIMUM 

POSITION 

DIFFERENCE 

(meters) 

HARR IS-PRI ESTER 

10.9 

108 

JACCHIA-ROBERTS 

11.2 

114 

MSIS 

11.0 

112 


NOTE: DURING THIS PERIOD, F^q.7 VARIED FROM 196.0 TO 
237.7 X 10“22 WATTS meter- 2 HERTZ“\ AND Kp 
VARIED FROM 1 TO 7. 
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SECTION 4 - CONCLUSION 


The results presented in this paper support the conclusion 
that, for satellites above 300 kilometers, the Harris- 
Priester, Jacchia-Roberts , and MSIS atmospheric density 
models all produce roughly similar density profiles and es- 
sentially comparable orbit determination results when the 
drag variation parameter is solved for and orbit quality is 
measured by adjacent arc overlap comparisons. It is impos- 
sible to predict which of the three models will produce the 
best fit or best predictions for any given orbit determina- 
tion arc. However, for some problem arcs, switching atmos- 
pheric models may result in marked solution improvements. 
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A GENERAL METHOD FOR COMPUTING THE TOTAL SOLAR 
RADIATION FORCE ON COMPLEX SPACECRAFT STRUCTURES 


F. K. Chan 

Scientific Analysts and Consiiltants, Inc. 
itlUi Heathfield Road, Rockville, Md, 20853 


ABSTRACT 

A general approach has been developed for computing 
the force due to solar radiation on an object of arbitrary shape. 
This method circumvents many of the existing difficulties in 
computational logic presently encountered in the direct analytical 
or numerical evaluation of the appropriate surface integral. It may 
be applied to complex spacecraft structures for computing the total 
force arising from either specular or diffuse reflection or even 
from non-Lambertian reflection and re-radiation. 
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SECTION 1 - INTRODUCTION 


The problem of computing the total force or total torque 
on a spacecraft due to solar radiation is, in general, very difficult. 
Mathematically, it requires the evaluation of a surface integral 
over only the illminated region of the surface. Even if the illu- 
minated region is known by some other means, the evaluation of 
the surface integral can still be very difficult analytically 
in the case of complex spacecraft structures. Moreover, if the 
illuminated region is not known a priori, the difficulties are 
compounded by having to determine self— shadowing. For non— convex 
objects, it is not trivially governed by a condition such as 
cos 0 0 where Q is the angle between the sun vector and the out- 

ward vector normal to the surface. In fact, the logic in the present 
methods becomes extremely complicated and is also not fool-proof. 
Additional difficulties are introduced by choosing a set of points 
(vertices) on the surface to form a network in approximating it; 
this inadvertently leads to book-keeping problems associated with 
selecting appropriate sets of points for computing surface elements. 

This paper presents a general method for performing the 
computations without encountering the difficulties described above. 

It does not attempt to evaluate the surface integral directly as 
it presents itself as done in the usual methods, but considers 
the same problem from a slightly different point of view which 
leads to the same results. 
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For convenience, let us use the following notation; 
u = unit vector along a specified direction 
V « any unit vector orthogonal to u, i.e. , u • v ■ 0 
w “ third unit vector forming orthogonal triad, i.e., w « u x 'v 
0 = origin of coordinate system 
P = any point on object's surface 
"t “ vector from 0 to P 
P'*= projection of P onto (v,w)-plane 
■r'= projection of "r onto (v,w)-plane 

(x,y,z) * reference orthogonal system for describing object's 
surface. 

In the present anailysis, it is advantageous to choose u to be opposite 

in direction to the incident solar radiation. (Alternatively, it can 

also be chosen to be in the same direction.) 

The vectors u and "r are known in the (x,y,z) system. 

In general, if V is any vector, then it may be more explicitly 

written in the (x,y,z)-space as V, v and has components V , V , V . 

^ * (x,y,z) X y z 

That is, we implicitly mean 






O-i) 


In view of the definition of the vector v, we may choose 

V = 0. Then, it may be shown that the other two components are given by 
z 

r/.a) 


=t ± 
r 








(/• 3 ) 
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From the definition of w, we obtain 


A 

W ss 




0'^) 


Therefore, any vector V, . can be transformed to V- . 

^x,y,z; ^u,v,w^ 

by the equation 




(t'S) 


where the transformation matrix T is given by 





<tT^ 





Then, using equation (1.6), the vector v, ^ is transformed 

\X,y, z; 

and we obtain 


> 

It 




(.t'f) 

/tv- == 


"f 


(h9) 

A^^J■ ~ 




r/<?) 


Consequently, the projection vector r' is simply given by 

^ t ^ ^ 

A »* t’” T^- A^ W 


(hto) 


The component r^ of the vector r is particularly important because, 
for a complex spacecraft structure, it can be used to yield the sxirface 
element which is directly exposed to solar radiation. This can be seen 
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as follows; For any given point on the (v,w)-plane (i.e., for any 

given vector "r'), the point on the spacecraft which is not shadowed 

is the one which has the maximum value of r^, independent of where 

the origin of the (u,v,w) coordinate system is chosen. (It would be 

the minimum value of r if the vector u had been chosen to be in 

u 

the same direction as the incident solar radiation.) To find the 
illuminated surface of the spacecraft, we proceed by dividing 
the (v,w)-plane into cells of area /iv Aw with cell centers (v,,w.). 

X J 

At these cell centers, the illuminated surface element is the one 
which has the maximum value of r^. In this way, the logic of deter- 
mining self-shadowing is extremely simple as compared to other methods 
which encounter considerable difficulty conceptually and computationally. 
Thus, given a vector "r' = (0,v^,w^), the vector r “ ^ ^ 

corresponding to the illuminated point is determined. It is then 
transformed to the (x,y,z) -space by the equation 

.T 


A 




= T A, 


•) 


( A /!) 


At this point r, the unit vector n, ^ normal to the surface 

(x,y,z)» (x,y,z) 

is then obtained by 

(hiz) 


A 

n 


|Vf| 

where ^(x,y,z) = 0 denotes the equation of the surface in a region 


A . 


containing r. For convenience, the direction of n is chosen such that 

O'/i) 


A A V 

n * u. ^ Q 


This choice of direction automatically makes n the outward unit normal 
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If the sxirface element belongs to a closed surface. Moreover, 
it establishes a direction for n in the case of a surface for which 
an outward unit normal is meaningless (such as a finite planar surface). 
The vector n^^ y. is then transformed to the (u,v,w)-space using 
the equation 


/\ 

•n 

(M.v; w) 


7 


A 


( A /-f } 


The cell (v^,w^) whose area is corresponds to 

a surface element whose area is denoted by AA. It is evident that 
we have 


AA 


AtrA ux 

u ) 


Air Aw 



(h/S) 


Therefore, the force AF exerted on this surface element is given by 


AF ■=> p AA 

where "p is the solar radiation pressure vector acting on the surface 
element. Under very general conditions of surface reflection and 
re-radiation, it can be shown that this pressure vector has the form 

f - - [c,S + f C,)ii] 0~n) 


where S is the solar radiation flux per unit area normal to the flux, 
c is the velocity of light, and ^is the angle between the sun vector 
and the normal to the surface element, i.e., 

CO & ^ n •li C(' 1^) 


The coefficients , 


and may change with time due to aging of 
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of the surface material by some complex process. 

For the case of specular reflection and diffuse (Lambertian) 


reflection, the Cj^’s are given by 





(/' 2 o) 

C - 'k 

0 - 2 ./) 


where »■ the fraction of incident radiation reflected specularly 
k^ = the fraction of incident radiation reflected diffusely 


by a Lambertian surface. 

It is to be noted that in equations (1.19) - (1*21 ), it is not 
implicitly assumed that the surface is radiating the entire energy 
incident on it, i.e., it is not necessary that we require the 
condition k^ + k2 ■ 1 in order to obtain these equations. 

For the case of speculen* reflection and non-Lambertian 
reflection and re -radiation a little consideration will reveal 

that the Cj^'s are given by 

= (I- fir) 

^ :i fiT 


c, = [r{/-^) + 




0-r)] 




Bf. ( V- ) 

where Y - the fraction of incident radiation reflected (specularly 
and otherwise ) 

^ =» the fraction of reflected radiation that is specular 
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Ej^ ■ non-Lambertian coefficients for front and back 
surfaces respectively 

e^f e^ • emission values for front and back surfaces 
respectively. 

In passing, it may be noted that we have the relations 

Moreover, it may be remarked that the form of equation (1.17) is 
valid for the more general non-Lambertian reflection and re-radiation 
which have a period of 7T in the azimuthal variable. In other words, 
Lambertian reflection means that the intensity I of the reflection 
is given by 

I I, 0'^7) 

Then, the case of non-Lambertian reflection and re-radiation expressed 
by equation (1.2U) would correspond to an intensity which is indepen- 
dent of the azimuthal variable <j> and is of the form 

1 - It f (^) (/>28) 

where implicitly we exclude the case of Lambert's law, i.e.. 

The even more general case means that we can have reflection and 
re-radiation for which the intensity is of the form 

I = 0‘io) 
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where 


’f (9j j>) - O'^O 

Finally, to compute the total force F due to solar radiation, 
we obtain from equations (1.15) - (1.19) the following expressions 





/’a ) 

( 0'-33) 

» - Av Auj ^ 21 [ 

C^) ^ J 

’ Ch34-) 

It is also trivial to compute the total torque M on the spacecraft 
by using the equations 

— ^ A i 

= /I X A F 


(/•3S-) 

M = 2^ 

but this will not be done here. 

} 

0>3i) 
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SECTION 3 - DISCUSSION 


It is obvious that the method just discussed does not 
encounter logic problems in determining self-shadowing. Moreover, 
because the points are first chosen on the projection plane, 

it circumvents the difficulties in book-keeping experienced in 
the other method of choosing vertices on the surface of the object. 
Furthermore, it does not require excessive core for storing the 
vertex data such as coordinates, area of surface element, noi*mal 
vector, solar incidence angle, etc. This advantage becomes evident 
by evaluating the expressions in equations (1.32) - (1.3U) using 
three accumulators (one for each force component), not having to 
store the set of points ^v^,w^| . Finally, if greater accuracy is 
desired, it suffices only to choose smaller values Av^w*, multiply 
the previous result by the factor Aw j * perform 

computations only for the additional points newly introduced into 
the set j^v^,Wj|. This advantage cannot be realized in the other 
method of choosing vertices on the surface of the object. In that 
case, in going to a refined model with additional vertices, it is 
necessary to perform the entire computations starting from the 
beginning each time. 
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SECTION It ~ CONCLUSION 

From the foregoing discussion, it may be concluded that 
the present method has the following advantages: 

1. It does not experience logic problems in determining 
self-shadowing. 

2 . It does not encounter the book-keeping problems arising 
in the case of choosing vertices on the svirface of the 
object. 

3 . It does not require excessive core for storing vertex 
data. 

ii. It can utilize previously obtained results in going to 
progressively more refined models. 
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SOLAR RADIATION FORCE MODELING FOR TORS ORBIT 

DETERMINATION * 

Taesul Lee, Michael J, Lucas, and Robert E. Shanklin, Jr. 

Computer Sciences Corporation 

ABSTRACT 

The relative orbit determination accuracies resulting from 
several TDRS models used for solar radiation force calcula- 
tions are evaluated. These models include spherical, single- 
plate, and restricted two-plate models. The plate models 
can be adjusted in both area and reflectivity through dif- 
ferential correction. The restricted two-plate model has an 
Earth-pointing plate and a solar plate; the orientation of 
the solar plate is restricted to rotation about an axis per- 
pendicular to the satellite's orbital plane. 

Simulated TDRS observations are generated from an ephemeris 
obtained using a 69-component TDRS model. These observa- 
tions are processed by least squares differential correction 
in order to find optimized parameters for the spherical, 
single-plate, and multi-plate models. The solutions for the 
parameters and the state vector are then used to generate 
ephemerides that are compared with the 69-component ephem- 
eris to estimate the expected orbit determination accuracies 
achievable with the various TDRS models. 


*This work was supported by the Operations Analysis Section, 
Operational Orbit Support Branch, Goddard Space Flight 
Center, National Aeronautics and Space Administration, under 
Contract NAS 5-24300. 
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SECTION 1 - INTRODUCTION 

A study of the solar radiation pressure (SRP) effect on 
orbit determination for a Tracking and Data Relay Satellite 
(TORS) has been carried out using simulated data. The TORS 
System consists of three geosynchronous satellites--TDRS 
East, TDRS West, and TDRS Spare--and one common ground 
tracking facility. These satellites will be placed in 
circular, nearly equatorial orbits at a height cf 36,000 kil- 
ometers above the surface of the Earth. The study is de- 
signed to determine whether a complex SRP model for a TDRS 
can be satisfactorily replaced by a simpler SRP model, such 
as a constant-effective-area model or a two-plate model. In 
addition, different tracking station configurations are used 
to investigate the possible dependence of the results on the 
tracking station geometry. 

A similar study carried out by Chan et al. (Reference 1) 
used a 69-Component TDRS SRP model and a two-plate model 
with four adjustable parameters. The adjustable parameters 
were determined by using a least squares procedure to mini- 
mize the position differences between two ephemerides, one 
obtained using the 69-component model and one obtained using 
the two-plate model. 

Another investigation related to the present study was 
carried out by Shanklin et al. (Reference 2) in which a 
constant-effective-area SRP model and a two-plate model were 
compared using real ATS-6 S-Band tracking data. This study, 
however, was somewhat incomplete due to the limited avail- 
ability of ATS-6 tracking data. The current study is an 
extension of that work and follows the same approach as that 
used in Reference 1 in constructing the TDRS SRP models. 

The current study, however, uses simulated bilateration and 
S-Band tracking data in the differential correction process 
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instead of position differences between the two ephemer ides 
as used by Chan. 

The study plan is as follows. First, a 69-component SRP 
model of a TORS, which is available in the Research and De- 
velopment version of the Goddard Trajectory Determination 
System (RDGTDS) , is used to compute a truth ephemeris, which 
is subsequently used to generate various types of simulated 
observations using the Mission Data Generation System 
(MDGS) . The MDGS produces raw data in a 75-byte format, and 
the Generalized Data Handler (GDH) converts these raw data 
into the 60-byte format for the Goddard Trajectory Determi- 
nation System (GTDS) . Second, these simulated data are used 
in regular GTDS Differential Correction (DC) Program runs to 
find optimized SRP parameters for the constant-effective- 
area model and for the two-plate model. The constant- 
effect ive-area model contains one adjustable parameter, and 
the two-plate model contains four adjustable parameters. 

Any combination of the four parameters of the two-plate 
model can be solved for in a given DC Program run. Third, 
ephemerides are generated using the final elements and SRP 
parameters obtained from the DC Program runs, and these eph- 
emerides are then compared with the original truth ephemeris. 

Brief descriptions of the TDRS solar radiation pressure 
models are given in Section 2 and generation of the simu- 
lated data is discussed in Section 3. The results of var- 
ious DC Program runs and ephemeris comparisons are presented 
in Section 4, and the conclusions are summarized in Sec- 
tion 5. 
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SECTION 2 - DESCRIPTION OF MODELS 

The 69-component model is composed of 69 distinctive parts. 
The components with relatively large areas are the two solar 
panels, whose normals make minimum angles with the satellite- 
Sun line; the antennas; the antenna feeds; and the top, bot- 
tom, and six sides of the main body (see Reference 1 for 
details) . 

The simplest SRP model used to approximate the 69-component 
model is the constant-effective-area model. In this model, 
the area for the SRP calculation is assumed to be constant 
and always normal to the satelli te-Sun line. The force due 
to the solar radiation pressure (Reference 3) is given by 



where v = eclipse factor 
a = constant area 

p = solar radiation pressure on a perfectly absorb- 
ing surface at the position of the satellite 

r\ = surface reflectivity 

~^Sun " unit vector along the satell ite-Sun line 

The solar radiation pressure is inversely proportional to 
the square of the distance from the Sun, and the eclipse 
factor, V, equals zero if the satellite is in the Earth's 
shadow and equals one if it is not. The right-hand side of 
Equation (2-1) represents the sum of two parts: the part 

due to the absorption of the solar radiation, which is pro- 
portional to (1 - r)) , snd the part due to the reflected 
radiation, which is proportional to 2n, This model is 
currently available in GTDS . 

The second model used to approximate the 69-component model 
is a two-plate model, which has an Earth-pointing plate and 
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a solar plate. The solar plate is hinged along an axis 
normal to the satellite's orbital plane and is always ro- 
tated about that axis so as to maximize the amount of sun- 
light falling on the plate. The force due to the solar 
radiation pressure for the two-plate model is given by the 
sum of four terms; 



where a = reference area 

P = solar radiation pressure on a perfectly absorb- 
ing surface at the position of the satellite 

ag = scale factor for the area of the Earth-pointing 
plate 

as = scale factor for the area of the solar plate 

riE = reflectivity of the Earth-pointing plate 

Hs = reflectivity of the solar plate 

^Sun = unit vector along the satell ite-Sun line 

R = unit position vector of the satellite 

Ng = unit vector normal to the sunny side of the 
solar plate 

In Equation (2-2) , the first term is due to the reflection 
by the Earth-pointing plate, the second term is due to the 
absorption by the Earth-pointing plate, the third term is 
due to the reflection by the solar plate, and the fourth 
term is due to the absorption by the solar plate. The two 
area scale factors, and a^, and the two reflectivities, 
hg and rigf are adjustable parameters. In a given DC Program 
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run, any combination or all of these four parameters can be 
solved for. Instead of a„, a , and an alternative 

£j hi S S 

set of four parameters, C3^.and may also be de- 

fined (and solved for) : 


^2 ~ 


C4 = a^d - 
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SECTION 3 - GENERATION OF SIMULATED DATA 

Tracking data for this study were generated using a satel- 
lite ephemeris tape obtained from a special RDGTDS Program 
load module that contains a 69-component TDRS model for SRP 
evaluation. This ephemeris tape was used by the MDGS Pro- 
gram to generate a second tape of raw range and Doppler sim- 
ulated data. This simulated data tape was used by the GDH 
Program to generate tracking data in a format appropriate 
for use in the GTDS two-plate load module. Two types of 
tracking data were generated in this manner: Applications 

Technology Satellite Ranging (ATSR) bilateration data and 
Unified S-Band (USB) two-way data. 

3 . 1 ATSR BILATERATION DATA 

ATSR bilateration data were generated using the ground 
station at White Sands, New Mexico, as the ATSR tracker and 
the ground stations at Mojave, California; Rosman, North 
Carolina; Madrid, Spain; Quito, Ecuador; and Santiago, 

Chile, as the ATSR ground transponders. Figure 1 shows the 
positions of these six sites in relation to the expected sub- 
satellite point for the relay satellite. 

Using these five tracker/ground transponder pairs, tracking 
data with the following characteristics were produced: 

o Frequency: 5600 MHz (C-Band) 

• Primary frequency offset: 5.8875 MHz 

« Transponder delay: 0 . 0 km 

» Tracking mode: satellite-to-ground phase-locked 

transponder 

® Major range tone/minor range tone: 100 kHz/8 Hz , 
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FIGURE 1. SITE LOCATIONS AND EXPECTED SUB^ATELLITE POINT 
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• uplink pilot frequency/downlink pilot frequency: 
,6150 MHz/4150 MHz 

• Doppler count mode: nondestruct 

Data were produced at a rate of six observations per minute 
for the first 25 minutes of each hour, starting at 0.0 hours 
on October 2, 1980, and ending at 0.0 hours on October 3, 
1980. Each tracker/ground transponder pair was enabled for 
tracking over the discrete time interval shown in Table 1. 

No observation corrections were applied and no observation 
noise was applied. 

3 . 2 USB TWO-WAY DATA 

USB two-way data (for which the receiving and transmitting 
sites are the same) were generated using the ground stations 
at Mojave, Rosman, Madrid, Quito, and Santiago. Tracking 
data with the following characteristics were produced: 

• Transmit frequency: 2106 MHz 

• Transponder delay: 0.0 km 

• Ranging equipment: Spaceflight Tracking and Data 

Network (STDN) Ranging Equipment (SRE) 

• Major range tone: 20 kHz 

Data were produced for the first 25 minutes of each hour, 
over the same time period, at the same rate, and with the 
same corrections that were used for the ATSR bilateration 
data. Each ground station was enabled for tracking over the 
discrete time interval shown in Table 2. 


3-9 



TABLE 1. TRACKING INTERVALS FOR ATSR TRACKER/ 
GROUND TRANSPONDER PAIRS 


TRACKER/GROUND 
TRANSPONDER PAIR 

MINUTES OF THE HOUR DURING 
WHICH THE PAIR IS ENABLED 

WHITE SANDS/ROSMAN 

00 TO 05 

WHITE SANDS/MOJAVE 

05 TO 10 

WHITE SANDS/QUITO 

10 TO 15 

WHITE SANDS/MADRID 

15 TO 20 

WHITE SANDS/SANTIAGO 

20 TO 25 


TABLE 2. TRACKING INTERVALS FOR USB GROUND STATIONS 


GROUND 

STATION 

MINUTES OF THE HOUR DURING 
WHICH GROUND STATION IS ENABLED 

ROSMAN 

00 TO 05 

MOJAVE 

05 TO 10 

QUITO 

10 TO 15 

MADRID 

15 TO 20 

SANTIAGO 

20 TO 25 
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SECTION 4 - DIFFERENTIAL CORRECTION SOLUTIONS AND 
~ EPHEMERIS COMPARISON RESULTS 

Differential correction solutions were obtained using dif- 
ferent SRP models, different types of simulated observa- 
tions, and different tracking station configurations. 

4.1 RESULTS OBTAINED USING BILATERATION DATA AND TWO GROUND 
TRANSPONDERS ^ ^ ^ ~ ^ ~ 

The results of DC Program solutions obtained using bilatera- 
tion range and Doppler data and five different combinations 
of solve-for parameters in the two-plate model are presented 
in Tables 3 and 4. The simulated bilateration data used 
were obtained using the TDRSS ground station at White Sands 
and two ground transponders at Rosman, North Carolina, and 
Mojave, California. The five different SRP options used were 

• Constant-effective-area model with C solved for 

R 

• Two-plate model with a_ and a solved for 

E s 

• Two-plate model with and ^2 solved for 

• Two-plate model with ?4 solved for 

• Two-plate model with ^ 2 > ^nd solved for 

The third option, in which and ^2 solved for, is 
equivalent to solving for and n^/ the scale factor and 
reflectivity of the Earth-pointing plate, respectively. 
Similarly, the fourth option is equivalent to solving for 
and n^. In this particular set of DC Program runs, the 
values of the SRP parameters in the two-plate model that 
were not solved for were set equal to zero. Thus, the third 
and fourth options discussed above actually represent single- 
plate models rather than two-plate models. 

An identical set of a priori elements, obtained from the 
truth ephemeris of the 69-component SRP model, was used for 
all of the options. It is seen from Tables 3 and 4 that the 
option of using the Earth-pointing plate alone gives the 
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TABLE 3. DIFFERENCES BETWEEN THE FINAL AND A PRIORI ELEMENTS 
(FINAL MINUS A PRIORI) 


CHANGES 

IN 

ELEMENTS 

RUN CONFIGURATION 

CONSTANT 
EFFECTIVE AREA 
SOLVED FOR 

TWO-PLATE MODEL 
AND as 
SOLVED FOR 

TWO-PLATE MODEL 
£1 and £2 
SOLVED FOR 

TWO-PLATE MODEL 
£3 AND £4 
SOLVED FOR 

TWO-PLATE MODEL 

ti, £2. and £3 
SOLVED FOR 

AX (meters) 

-5.81 

-3.42 

-24.16 

-1.62 

-2.63 

AY (meters) 

12.62 

-9.57 

- 18.43 

-8.59 

-7.61 

AZ (meters) 

-24.80 

16.07 

-109.73 

-38.56 

11.17 

AX (cm/sec) 

0.052 

0.036 

0.250 

0.008 

0.029 

AY (cm /sec) 

-0.041 

-0.024 

-0.221 

-0.005 

-0.025 

AZ Icm/sec) 

-0.118 

-0.138 

0.243 

-0.141 

- 0.084 


NOTES: 1. THE SAME SET OF A PRIORI ELEMENTS WAS USED FOR ALL DC PROGRAM RUNS. 

2. THE QUANTITIES ttg AND «5 DENOTE SCALE FACTORS FOR THE AREAS OF THE EARTH-POINTING PLATE AND THE SOLAR 
PLATE, RESPECTIVELY. THE PARAMETERS J,, £2- ^3- AND £4 ARE DEFINED AS FOLLOWS: £2 = “g " “S’'S- 

£4 = «s <1 - ’IS>- WHERE i|E AND 7,3 DENOTE THE REFLECTIVITY OF THE EARTH POINTING PLATE AND THE SOLAR PLATE, 

RESPECTIVELY. 
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TABLE 4. DC PROGRAM STATISTICS AND SRP PARAMETERS SOLVED FOR 




RUN CONFIGURATION 


PARAMETERS 

CONSTANT 
EFFECTIVE AREA 
SOLVED FOR 

TWO-PLATE 
MODEL 
Off AND urs 
SOLVED FOR 

TWO-PLATE 
MODEL 
AND h 
SOLVED FOR 

TWO-PLATE 
MODEL 
£3 AND ^4 
SOLVED FOR 

TWO-PLATE 

MODEL 

and |3 

SOLVED FOR 

WEIGHTED RMS 
STANDARD DEVIATION 

0.0558 

0.0346 

0.3238 

0.0546 

0.0329 

RANGE Imeters) 

0.514 

0.584 

2.206 

0.481 

0.593 

DOPPLER (millihertz) 

0.914 

0.473 

5.416 

0.898 

0.431 

SRP PARAMETERS 
SOLVED FOR 

Cr = 1.38 

= 0.281 
as = 1.219 

= 1.971 
= 0.551 

l3 - 19.482 

= -37.588 

= 0.175 
^2 = 0.146 
^3 = 0.602 
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poorest results, whereas the other options all give compar- 
able results. Similar conclusions are supported by Fig- 
ures 2 and 3, which represent 24-hour ephemeris comparison 
results between the original 69-component ephemeris and the 
ephemerides obtained using the DC Program solutions for dif- 
ferent SRP options. The results obtained using the second 
option, in which and were solved for, are not shown 
because they are very similar to the results obtained using 
the fourth option. Only the along-track and cross-track po- 
sition differences are shown in Figures 2 and 3 ,■ because the 
radial position differences were much smaller than the along- 
track or cross-track position differences. 

The single-plate option using the Earth-pointing plate alone 
gives the worst position errors. The single-plate option 
using the solar plate alone gives significantly better re- 
sults. in fact, the option using the solar plate alone 
gives the smallest along-track position differences of all 
the different options used. 

There are two features worth mentioning. First, there is no 
significant difference between the constant-effective-area 
model and the more complex two-plate model options. Second, 
in all cases studied, there are quite sizable cross-track 
position differences, equal to or larger than the along- 
track differences. 

in order to examine the influence of the tracking geometry 
on the orbit determination results, a different pair of 
ground transponders (Rosman and Santiago) was used for the 
same series of DC Program solutions discussed above. Ephem- 
eris comparison results obtained using these differential 
correction solutions were then compared with the correspond- 
ing results obtained using the pair of ground transponders 
at Rosman and Mojave; the only significant difference between 
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ALONG-TRACK DIFFERENCES (niutBisI 



: CONSTANT EFFECTIVE AREA VERSUS 69-COMPONENT 
: TWO-PLATE (THREE PARAMETERS SOLVED FOR) VERSUS 69-COMPONENT 
■«' ■ : SINGLE SOLAR PLATE VERSUS 69-COMPONENT 

; SINGLE EARTH-POINTING PLATE VERSUS 69-COMPONENT 

FIGURE 2. ALONG-TRACK POSITION DIFFERENCES FROM 24-HOUR EPHEMERIS 
COMPARISON RESULTS 
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TIME (hours) 

CONSTANT EFFECTIVE AREA VERSUS 69-COMPONENT 
TWO-PLATE (THREE PARAMETERS SOLVED FORI VERSUS 69-COMPONENT 
SINGLE SOLAR PLATE VERSUS 69-COMPONENT 
SINGLE EARTH-POINTING PLATS VERSUS 69-COMPONENT 

FIGURE 3. CROSS-TRACK POSITION DIFFERENCES FROM 24-HOUR EPHEMERIS 
COMPARISON RESULTS 





the two sets of results was in the cross-track position dif- 
ferences. The maximum cross-track position differences ob- 
tained using the Rosman and Santiago ground transponders 
were found to be less than 10 meters, whereas the corre- 
sponding differences obtained using the Rosman and Mojave 
ground transponders were larger than 20 meters. 

4.2 RESULTS OBTAINED USING S-BAND RANGE DATA AND TWO GROUND 

TRACKING STATIONS ^ ^ 

Differential correction solutions for a 24-hour TDRS arc 
were obtained using S-Band range data and two different 
tracking station configurations. In the first set of solu- 
tions, the two ground stations at Rosman and Mojave were 
used, and in the second set of solutions, the two stations 
at Rosman and Santiago were used. The results of 24-hour 
ephemeris comparisons are summarized in Figures 4 and 5. it 
is seen from Figures 4 and 5 that the results obtained using 
S-Band range data are generally worse than the corresponding 
results obtained using bilateration data. The along-track 
position differences shown in Figure 4 indicate that the 
Rosman/Mo jave configuration gives somewhat better results 
than does the Rosman/Santiago configuration. in the case of 
the cross-track position differences shown in Figure 5, the 
situation is reversed; the Rosman/Santiago configuration 
gives somewhat better results than does the Rosman/Mo jave 
configuration. 

4-3 RESULTS OBTA INED USING MORE THAN TWO GROUND TRACKING 

stations' ^ 

The same 24-hour TDRS arc studied in Sections 4.1 and 4.2 
was used in a set of DC Program runs using more than two 
ground tracking facilities. in the case of bilateration 
data, three ground transponders, located at Mojave, Santiago, 
and Madrid, and five ground transponders, located at Mojave, 
Santiago, Madrid, Rosman, and Quito, were used. Ephemeris 
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TRACKING CONFIGURATION 
ROS, GDS ROS, AGO 

— o — o— — •- 



CONSTANT EFFECTIVE-AREA vs 69-COMPONENT 
TWO-PLATE l {3 AND {4 SOLVED FORI vs 69-COMPONENT 
TWO-PLATE It^, ( 2 ' '3' ^4 SOLVED FORI vs 69-COMPONENT 


FIGURE 4. ALONG-TRACK POSITION DIFFERENCES FROM 24-HOUR EPHEMERIS 
COMPARISON RESULTS 


CROSS TRACK DIfFERENCES (mmets) 




CONSTANT-EFFECTIVE-AREA vs 69-COMPONENT 
TWO-PLATE (fj AND SOLVED FORI vs 69-COMPONENT 
TWO-PLATE (?,, E^. s'3 . and E 4 SOLVED FOR) vs 69-COMPONENT 



FIGURE 5. CROSS-TRACK POSITION DIFFERENCES FROM 24-HOUR EPHEMERIS 
COMPARISON RESULTS 
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comparison results obtained using the three ground trans- 
ponders were similar to the results obtained using the five 
ground transponders. Typical along-track, cross-track, and 
radial position differences were 6.0, 1.0, and 1.0 meters, 
respectively. No significant difference was found among the 
different models used for the solar radiation pressure com- 
putation as long as the initial state vector and the solar 
radiation pressure parameters were solved for. 

Similar analyses were carried out using more than two S-Band 
tracking stations. Two sets of differential correction so^ 
lutions were obtained using three tracking stations at 
Mojave, Madrid, and Santiago and four tracking stations at 
Mojave, Rosman, Madrid, and Santiago. Ephemeris comparison 
results obtained using these differential correction solu- 
tions are summarized in Tables 5 and 6 . There is no essen- 
tial difference between the results obtained using three 
tracking stations and the results obtained using four track- 
ing stations. These results show a significant improvement 
over the corresponding results obtained using only two 
S-Band tracking stations. Cross-track position differences 
were reduced by almost a factor of 10 and along-track dif- 
ferences were also substantially reduced. However, none of 
the results obtained using S-Band tracking data were as good 
as the corresponding results obtained using bilateration 
data . 
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TABLE 5. CROSS-TRACK AND ALONG-TRACK POSITION DIFFERENCES OBTAINED 
USING THREE USB GROUND STATIONS (MAD, AVE, AGO) 


EPHEMERIDES COMPARED 

MAXIMUM CROSS-TRACK 
DIFFERENCE (meters) 

MAXIMUM ALONG TRACK 
DIFFERENCE (meters) 

CONSTANT-EFFECTIVE-AREA vs 69-COMPONENT 

3.1 

24.6 

TWO-PLATE («£ AND ug SOLVED FOR) vs 69-COMPONENT 

2.8 

23.5 

SINGLE SOLAR PLATE (^3 AND ^ SOLVED FOR) vs 69-COMPONENT 

7.1 

22.2 

TWO PLATE (£ 1 , ^ 3 - AND £4 SOLVED FOR) vs 69-COMPONENT 

6.5 

20.2 
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TABLE 6. CROSS-TRACK AND ALONG-TRACK POSITION DIFFERENCES OBTAINED 
USING FOUR USB GROUND STATIONS (MAD, AVE. AGO, ROS) 


EPHEMERIDES COMPARED 

MAXIMUM CROSS-TRACK 
DIFFERENCE (meters) 

MAXIMUM ALONG-TRACK 
DIFFERENCE (meters) 

CONSTANT-EFFECTIVE AREA vs 69 COMPONENT 

3.3 

26:2 

TWO PLATE («£ AND ug SOLVED FOR) vs 69-COMPONENT 

3.2 

25.9 

SINGLE SOLAR PLATE (Jg AND ^ SOLVED FOR) vs 69-COMPONENT 

4.4 

26.7 

TWO PLATE (1^, ^ 2 - ^ 3 - ^ND ^4 SOLVED FOR) vs 69 COMPONENT 

4.2 

25.0 
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SECTION 5 - CONCLUSIONS 

A study of solar radiation pressure (SRP) as it affects TDRS 
orbits was performed using simulated bilateration data, sim- 
ulated direct two-way data, and various ground station con- 
figurations. Orbit determination results obtained using 
constant-effective-area and two-plate SRP modeling were com- 
pared with each other and. with an ephemeris obtained using a 
69-component TDRS SRP model. The conclusion of this study 
can be summarized as follows: 

• The constant-effective-area solar radiation pres- 
sure model and the two-plate model give essentially 
the same quality results when both the state and 
the SRP parameters are solved for. The maximum 
position differences between the 69-component model 
truth ephemeris and an ephemeris determined using 
solved-for elements and SRP parameters can be re- 
duced to less than 10 meters if proper bilateration 
tracking configurations are used in solving for the 
elements and the SRP parameters. 

a When using only two ground tracking facilities, the 
Rosman/Santiago combination gives smaller cross- 
track position errors than does the Rosman/Mo jave 
combination . 

®, Results obtained using three ground tracking facil- 
ities (located in a triangular configuration) are 
significantly better than the corresponding results 
obtained using two ground tracking facilities. 

® Results obtained using more than three ground 

tracking facilities are of essentially the same 
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quality as the results obtained using three ground 
tracking facilities. 

• Bilateration data appear to give better orbit de- 
termination results than S-Band tracking data. 
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PRECISION ORBIT COMPUTATIONS FOR AN OPERATIONAL ENVIRONMENT 

C. E. Doll, Goddard Space Flight Center 
David F. Eggert, Computer Science Corporation 
Richard L. Smith, Computer Science Corporation 

ABSTRACT 

Analyses have been performed at the Goddard Space Flight Center (GSFC) to 
establish the operational procedures that would be required to provide pre- 
cision orbit computations to meet current and future operational requirements 
set forth by different NASA projects. Taking advantage of the improvements 
to the earth's gravitation field and tracking station coordinates, an orbital 
computational consistency of the order of 5 meters were achieved for total 
position differences between orbital solutions for the Seasat and GEOS-3. 

The main source of error in these solutions has been in the mathematical models 
that are required to generate these results, i.e., gravitation, atmospheric 
drag, etc. Different earth's gravitation fields and tracking coordinates have 
been analyzed and evaluated in obtaining these computational results. 

Comparisons and evaluations of the Seasat results have been obtained in terms 
of different solution types such as the Doppler only. Laser only, Doppler and 
Laser, etc. Other investigation using the Seasat data have been made in 
order to determine their effect on the computational results at this partic- 
ular level of consistency. 
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INTRODUCTION 

It is expected that in the next few years that NASA missions will require 
additional computational precision in determining spacecraft position in 
order to support both project and scientific requirements. In order for the 
Goddard Space Flight Center to support these NASA mission in a precision 
orbit computations environment both methods and techniques for computations 
and operational procedures must be established. 

The definitive orbit computations requirements for the Seasat mission were 
the most accurate in terms of consistency between orbital solutions that had 
been performed at the GSFC for any given mission prior to its launch in June 
1978 by the Operations Support Computing Division (OSCD). The computations 
requirements set forth by the Seasat Project was to maintain a maximum devi- 
ation of 65 meters between orbital solutions for the mission lifetime. With 
these project requirements, the OSCD established the computational techniques, 
the operational procedures and the tracking data distribution in order to ful- 
fill these commitments. 

Due to the amount and distribution of USB/SRE and Laser tracking data required 
to support definitive orbit computations and precision orbit computations for 
the Seasat mission, the OSCD has taken the initiative to determine what level 
of consistency between orbital solutions can be reached for an operational 
environment. The results of these investigations for the Seasat and GEOS-III 
missions are based on the mathematical models and station geodetics that have 
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Introduction (continued) 

been established at GSFC by the Geodynamics Branch, The computational pro- 
cedures and observational tracking data distributions have been established 
through the analyses which have been performed for each of the satellites. 

The information in this particular report is presented in three different 
areas, the method for precision orbit computations, Seasat precision com- 
putations and GEOS-III precision orbit computations. 



METHODS FOR PRECISION COMPUTATIONS 
Orbit Determination Procedure 

The computations of the precision orbits for both Seasat and GEOS-III were 
performed at the GSFC on the 360 computer complex using the Goddard Tra- 
jectory Determination System (GTDS). GTDS has the capability to perform 
orbit determinations and generate spacecraft ephemeris data in the form of 
position and velocity to different levels of consistency based on force 
model representations, station geodetics and tracking data distributions. 

The orbital solutions obtained for Seasat and GEOS-III from GTDS used 
Cowell's method of integration for the equations of motion and the vari- 
ational equations and a least squares adjustment technique for the improve- 
ment of orbital parameters. The earth's gravity field, the solar gravita- 
tional perturbations, the lunar gravitational perturbations and the solid 
earth tidal perturbations are modeled for these orbital computations. In 
addition. The nonconservative forces of solar radiation pressure and atmos- 
pheric drag have been modeled. It should be stated that the JPL planetary 
ephemeris DE-96 was adopted for these computations along with the BIH polar 
motion and the UTl and A. 1 corrections. 

The Seasat and GEOS-III spacecraft were modeled in the GTDS as specularly 
reflecting spheres. In the precision orbit computations for Seasat a drag 
coefficient for each data arc was solved for. 

In addition, an analysis was performed to determine the best integration step 
size for the equations of motion and the variational equations and in obtain- 
ing orbital solutions which are consistent in terms of numerical processes. 
The integraton step size which was established for Seasat and GEOS-III was 
45 seconds. 
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Physical Parameters, Environmental Parameters and Tracking Station Geodetics 
For Precision Orbit Determination 

In obtaining the orbital solutions for the Seasat and GEOS-III in the pre- 
cision orbit computations environment different sets of physical and environ- 
mental parameters and station geodetics were used and evaluated. One of the 
fundamental capabilities that exist in GTDS is its capability to make use of 
different size gravitational models along with other parameters, which is 
essential in an operational environment. In this investigation the three 
earth's gravitational fields which were used and evaluated were the GEM 9, 

GEM lOB, and the PGS 1040. These three gravitational fields were determined 
at the GSFC using observational tracking data from both NASA and non-NASA 
stations and global gravimetric data while making use of the research and 
development orbit computations system GEODYN. When a specific gravitational 
field is used for orbit computations then the earth's gravitational constant 
(GM), the mean equatorial radius of the earth (ag) and the earth's inverse 
flattening factor (1/f) must be properly specified. These particular parameters 
for each of the three gravitational fields are listed in Table 1. The orbital 
and physical parameters that were used in this investigation are listed in 
Table 2. It should be understood that in the computations for the noncon- 
servative forces of drag and solar radiation that both spacecrafts were assumed 
to have a spherical shape, although this is usually an extreme Idealization. 

Through the analysis and evaluations which have been performed in this invest- 
igation for precision orbit computations, it has become apparent that good 
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(Physical Parameters, etc. , continued) 


or precise station geodetlcs are very essential In obtaining specific levels 
of consistency between orbital solutions. The evaluations which have performed 
Indicates that the quality of station geodetlcs are not as Important at the 20 
to 40 meter level of consistency between orbital solutions as they are at the 
5 to 15 meter level of consistency between solutions. Therefore, the station 
geodetlcs which have been used for the precision orbit computations for both 
Seasat and GEOS-llI are the coordinates which have been derived by J. Marsh 
of the GSFC which are given In Table 3. It should be pointed out that 
selected code letters are assigned to specific stations In order to represent 
that station on the tracking data distribution figures that are presented In 
Figures 1 through 3. 
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SEASAT PRECISION ORBIT COMPUTATION 


Observational Tracking Data for Seasat 

The observational tracking data used for precision orbit computations for 
Seasat were a combination of USB/SRE range rate data from STDN and Laser data 
from STDN and SAO. The USB/SRE range rate data provided the strong global 
coverage both in terms of geographical distribution and in time. The Laser 
observational tracking data provided strength in terms of accuracy for the 
precision orbit computations. 

An analyses of both the USB/SRE range rate data and the Laser data in terms 
of distribution and time provided two specific time intervals, September 19 
through September 26, 1978 and August 8, 1978 through August 15, 1978 over 
which the precision orbit computations were performed. The amount of obser- 
vational tracking data during these two particular time intervals contained 
approximately 20 passes of USB/SRE data and 12 passes of Laser data for each 
typical twenty-four hour interval. Figures 1 and 2 give the station and 
data distribution for the September 1978 period and the August 1978 period. 

Orbital Analyses for Seasat 

In determining the consistency between orbital solutions to the 1 to 5 meter 
level for the Seasat spacecraft, a number of gravitational field models, 
station geodetics and integration step size were evaluated. Through these 
evaluations with the use of GTDS, it has been established that the PCS- 1040 
gravitational field and the station geodetics, which have been designated Marsh 
II, have given the best results in terms of consistency between orbital solu- 
tions. The PGS-1040 gravitational field and the Marsh II station geodetics 
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(Orbital Analysis for Seasat - continued) 

have been determined at GSFC through the use of GEODYN. It should be pointed 
out that in the determination of the PGS-1040 gravity field that both Laser 
and USB/SRE observational tracking data from the Seasat spacecraft were used. 

The length of the observational data arc was thirty hours for the orbital 
solutions which were determined for this investigation. In order to deter- 
mine the consistency between successive orbital solutions for the Seasat 
spacecraft a six-hour interval was established as the time frame over which 
the consistency was to be determined. The maximum difference in a given 
six-hour overlap interval between two successive orbital solutions in terms 
of spacecraft position is the measure of consistency which has been deter- 
mined by this process. 

The orbital solutions for the Seasat spacecraft using only the USB Doppler 
tracking and the additional techniques for computations in the September 
and August 1978 time frames are given in Tables 4 and 8. Information per- 
taining to the individual solutions are given in these tables including the 
rho one solve-for parameter, which is equivalent to a density correction for 
each of the Seasat orbital solutions. In addition, the maximum discon- 
tinuties between successive solutions for each specific six— hour overlap 
interval are presented in terms of radial, cross track and along track dif- 
ferences. The results of this analysis indicate that using the Doppler only 
that an average 10-meter level of consistency for the September 1978 time 
frame can be obtained while for the August 1978 time frame only a 13-meter 
level of consistency was obtained. These results indicate that the 5-meter 
level of consistency between the orbital solutions is difficult to obtain 
using only USB Doppler data. An assessment of these results would indicate 
that there should be no problem with the number of tracking passes in the 
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Orbital Analyses for Seasat (continued) 


individual solutions although the distribution of passes within the solutions 
could cause problems. It is felt that the mathematical modeling or the com- 
putational procedures should not cause problems in achieving the 5-meter level 
of consistency. 

The next set of orbital solutions for Seasat were computed based on Laser 
tracking data only and the results of these computations are given in Tables 
5 and 9. Information pertaining to these computations for the individual 
solutions are given in these tables including the rho one solve-for parameters 
The maximum discontinuities between successive orbital solutions for each 
specific six-hour overlap Interval are presented. The results of this analy- 
sis indicate that using the Laser tracking data by itself that an average 
4.4 meter level of consistency can be obtained for the September 1978 time 
frame while for the August 1978 time frame only an 8.8-meter level of 
consistency was obtained. These results indicate the 5-meter level of con- 
sistency between individual solutions can be obtained when using only Laser 
tracking data for certain time frames during the Seasat satellite lifetime. 
Again, an assessment of these results would Indicate that since the mathe- 
matical modeling and the computational procedures are the same then the 
differences in the August and September 1978 time frames has to be in an- 
other area. The only other area where differences can be attributed has to 
be in the Laser tracking data, in other words the distribution of the data 
or the quality of data. 
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Orbital Analyses for Seasat (continued) 


Another set of orbital solutions for Seasat were determined based on Laser 
and USB Doppler tracking data and the results of these computations are 
given in Tables 6 and 10. The information pertaining to these computations 
are given in these tables, including the rho one solve-for parameters. The 
maximum discontinuities between successive orbital solutions for each speci- 
fic six-hour overlap interval are also presented in these tables. The 
results of this analysis indicate that using both the. Laser and USB Doppler 
tracking data that an average 3.6-meter level of consistency was obtained 
for the September 1978 time frame while for the August 1978 time frame 
only a 7.4-meter level of consistency was obtained. These results indicate 
that making use of the combination of Laser and USB Doppler tracking data 
gives a little better overall consistency between successive solutions than 
when using the Laser observations only. Since the mathematical modeling 
and the computational procedures were the same then the slight improvements 
comes from the strength of more comprehensive distribution of observational 
tracking data throughout the individual orbital solutions. 

Further analysis was performed to determine the affect of having equal number 
of observations per pass for both the Laser and USB Doppler tracking data in 
determining each orbital solutions and the level of consistency for the 
September 1978 time frame. The results of these individual orbit computations 
are given in Tables 6 and 7 along with the rho one solve-for parameters. The 
maximum discontinuities between successive orbital solutions for each six-hour 
overlap interval are also presented in these tables. The results of this 
analysis indicate that making use of the observational tracking data in this 
manner and using the same mathematical modeling and computational procedures 
an average of 4. 1 meter level of consistency was obtained. This result of 
4.1-meter level of consistency obtained in this process and the other average 
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Orbital Analyses for Seasat (continued) 


values of 3.7- and 4.4-meter levels of consistency obtained when using Laser 
and USB Doppler data in another process of observations selection and using 
Laser data by itself are basically the same. In other words, at this 
particular level of consistency it is difficult to indicate in terms of an 
average value, which are the better results. 
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GEOS-III PRECISION ORBIT COMPUTATIONS 
Observational Tracking Data for GEOS-III 

GEOS-III orbital solutions were calculated for a period extending from 
February 23, 1976, to March 2, 1976. The available unified S-band range 
and range-rate data is shown in Figure 3. Only the range-rate data were used 
for the solutions described here. Unlike the tracking data distribution 
for Seasat, the GEOS-III tracking data distribution is not uniform, having 
intense tracking about once a day, and very little tracking at other times. 

On the average, there is available slightly less than one pass of tracking 
per orbital revolution. 

Orbital Analysis for GEOS-III 

Orbital solutions for GEOS-III were calculated using GTDS and the Goddard 
Earth Model lOB (GEMIOB) gravity model. This gravity model is based, in part, 
on GEOS-3 altimetry data. Since the altitude of GEOS-III is about 50 kilo- 
meters greater than that of Seasat, the orbital effects of atmosphere drag 
are significantly smaller. Unlike Seasat, estimation of the drag parameter 
does not sppear to affect the accuracy of differential correction solutions. 
The GEOS-III solutions were calculated by solving only for the spacecraft 
state vector at epoch. 

The GEOS-III solutions were 30 hours in length, each solution overlapping 
neighboring solutions by six hours. Because ephemeris comparisons In the 
solution overlap Intervals are used for orbital accuracy estimates and because 
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Orbital Analysis for GEOS-III (continued) 


of the strongly periodic characteristic of the tracking schedule, It might be 
expected that the overlap comparisons could be affected by the placement of 
the overlap Interval relative to the periods of Intense tracking. If the 
overlap intervals coincided with the intense tracking periods it might be 
expected that the ephemeris differences would be lower than If the overlap 
intervals were located in periods of little tracking. 

In order to examine this possible effect, the solution intervals were placed 
in time two different ways. In the first scheme, the epochs of each 30-hour 
solution were located at 15^ on successive days. This procedure puts the 
periods of intense tracking into the six-hour solution overlap intervals, 
and each soluton has strong tracking at its start and end, but little in 
between. The second scheme placed the epochs at on successive days. This 
placed the intense tracking in the middle of each solution, with very little 
in the overlap intervals. 

GEOS-III orbital solutions, along with the ephemeris overlap comparisons that 
were calculated using these two approaches are summarized in Tables 11 and 12. 
In these tables, the tracking observations for each solution are separated 
into two categories (indicated by the diagonal line) because of slightly 
different tracker types; this is not relevant for this study. The orbital 
fits, as indicated by the weighted RMS, (the assigned range-rate standard 
deviation was 2.0 centimeters per second) were about the same, overall, for 
the 0^ and 15^ solutions. Similarly, the standard deviations of the solution 
residuals were about one centimeter per second for each set of solutions. 
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Orbital Analysis for GEOS-III - continued 


The ephemeris overlap differences for both sets of solutions are also quite 
similar. The maximum total differences average about 7 meters for both the 
0 ^ and 15^ solutions. Also the maximum cross-track differences average about 
6 meters for both sets of solutions. On the other hand, the radial and along- 
track differences for the two sets of solutions are distinct. For the 15^ 
solutions, the maximum radial differences and the maximum along-track differ- 
ences average to 0,5 and 2.4 meters, respectively. For the 0^ solutions, the 
corresponding averages are 1.0 and 4.9 meters. Thus, the placement of the 
intense tracking at the end of the solution intervals, rather than the middles, 
reduced the along-track and radial differences by about a factor of two. 

This reduction in along-track and radial differences, and presumably, a 
corresponding reduction in along-track and radial orbit error may be explained 
as follows. It is well known that radial and along-track orbit displacements 
are coupled together in the equations of motion; thus it is natural that 
changes in along-track and radial orbit error should be correlated. Placement 
of the intense tracking at the ends of a solution interval causes the orbit 
solution to better average out along-track and radial force modeling errors, 
leading to smaller peak radial and along-track orbit errors than if the 
tracking data was concentrated in the middle of each solution, leaving both 
ends of a solution "floating". 
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COMPARISONS OF VARIOUS SETS OF TRACKING STATION COORDINATES 

The GEOS-III solutions described in the previous section were calculated 
using tracking station coordinates derived by J. Marsh of GSFC. Corres- 
ponding GEOS-III orbital solutions were calculated using three other sets 
of tracking station coordinates. These three sets are NASA Spacecraft 
Tracking and Data Network coordinates (STDN), GEM9 coordinates, and World 
Geodetic System (Geoceiver) WGS(G) coordinates. 

The STDN coordinates are those used for GSFC operational orbit determination 
(Reference A). The GEM9 coordinates were derived as a part of the GEM9 
and GEMIO gravity models (Reference B). The WGS(G) coordinates for the NASA 
S-band tracking stations were specially derived for this study. These 
station coordinates were based upon coordinates of nearby geoceivers. 

GEOS-III orbital solutions using the STDN, GEM9, and WGS(G) station coordin- 
ates are summarized in Tables 13, 14, and 15 respectively. These solutions 
were calculated using the same GTDS input parameters, except for station 
coordinates as the solutions in Table B (15^ epochs). Thus, comparisons 
among the results in these four tables are a direct comparison of the effect 
of various sets of tracking station coordinates. (The value of the semi- 
major axis of the earth, used for evaluation of the gravity force was 
slightly different for the solutions calculated using Marsh coordinates. 
Subsequently, tests shovred the effect of this change negligible for these 
comparisons.) 
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Comparisons of Various Sets of Tracking Station Coordi nates (continued) 

None of the three additional sets of station coordinates performed as well in 
these solutions as the Marsh coordinates. In the order of increasing weighted 
RMS residuals and increasing overlap differences, these three sets of coor- 
dinates are ordered as follows: WGS(G), GEM9, and STDN. In the case of the 

STDN coordinates, the maximum radial differences average to 4.2 meters, while 
the total differences average to 21 meters. These results are consistent 
with the position differences of the GEOS-III tracking stations in the Marsh 
and STDN coordinates, which are typically 15 to 25 meters. 
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CONCLUSIONS 


The results of this study have shown that orbital consistency at the five- 
meter level can be obtained for Seasat and GEOS-III using the operational 
Goddard Trajectory Determination System. The attainment of this orbital 
consistency level requires the use of the most precise gravity models and 
tracking station coordinates that are currently available. For Seasat, 
the use of Laser range tracking data vras found to increase the level of 
orbital consistency when used alone or in combination with the unified S- 
band range-rate tracking data. For GEOS-III, the use of the unified S-band 
tracking data alone produced orbital consistency of the order of five 
meters. 
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Table 1 


Physical, Geophysical, and Astronomical 
Parameters Used 


QUANTITY 

VALUE 

UNIVERSAL CONSTANT 0f= GRAVITATION (G) 

6.673 X 10“^^ KM^ S“^KG~’ 

ASTRONOMICAL UNIT 

1.495978930 x 10® KM 

SOLAR MOMENTUM FLUX DENSITY 

4.5 N KM~^ 

EARTH GRAVITATIONAL CONSTANT (GM) 

3.9860064 x 10® KM® S~® (GEM 9) 
3.9860064 x 10® KM® (GEM10B) 

3.9860062 x 10® KM® S~® (PGS 1040) 

EARTH MEAN EQUATORIAL RADIUS (a ) 

e 

6378.140 KM (GEM 9) 

6378.139 KM (GEM10B) 

6378.140 KM (PGS 1040) 

EARTH INVERSE FLATTENING FACTOR (1/f) 

298.250 (GEM 9) 
298.257 (GEM 108) 
298.257 (PGS 1040) 

SPEED OF LIGHT (c) 

2.997925 X 10® KM S~’ 
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Table 3. Marsh II Tracking Station Coordinates 


STATION. 

GEODETIC 

LATITUDE" 

GEODETIC 

LONGITUDE 

■ HEIGHT ABOVE 
SPHEROID 
(m) 

CODE 

ACN3 

-7°57'17''.289 

345°40'22".186 

534.33 

A 

AG03 

-33°09'03",946 

289°20'00".558 

717.59 

B 

8DA3 

32°2r04".533 

295°20'31".325 

-30.10 

C 

ETCA 

38°59'54".171 

283°09'28".749 

12.35 

D 

GDS3 

35°20'31".789 

243°07'35''.311 

919,69 

G 

GDS8 

3S°20'29".495 

243°07'34".792 

925.69 

H 

GWM3 

13°18‘38".243 

144°44'12".465 

133.05 

1 

HAW3 

22°07'34".681 

200‘’20'05".231 

1148.56 

J 

MAD8 

40°27'19”.5o3 

355°49'53".216 

819.66 

K 

MIL3 

28°30'29".250 

279°18'23".625 

-38.24 

L 

ORR3 

-35°37'40".410 

148°57'25".169 

934.39 

N 

QUIS 

-Q°37'18",967 

281°25'10".404 

3578.86 

0 

ULA3 

64°58'19".233 

212°29'13'',235 

333.90 

Q 

MAD3 

40°27'22".248 

355°49'49".163 

816.30 

R 

MILA 

28°30'29".318 

279°18'25".474 

-42.40 

S 

AREL 

-16°27'56".70S 

288°30'24".533 

2475.99 

a 

BDAL 

-32°21'13".767 

295°20'37”.890 

-36.87 

b 

GTKL 

21°27'37".770 

288°52'04",972 

-32.36 

c 

HOPL 

31°41'03".201 

249°07'18".798 

2334.76 

d 

KOOL 

52®10'42".215 

5°48'35".055 

75.0 

e 

NATL 

-5°55'40'M45 

324°50'07".165 

22.70 

f 

ORRL 

-35°37'29".741 

148°57'17".133 

932.45 

g 

RAML 

28°13'40".630 

279°23'39".244 

-37.24 

h 

SNDL 

32°36'02".628 

243°09'32".737 

975.00 

i 

STAL 

39°01'13".359 

283°10'19".751 

47.00 

i 


Preference spheroid; semimajor axis, 6378.155 km. inverse flattening factor, 298.255. 
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TABLE 4 


SATELLITE AND TIME PERIOD SEASAT - September 1978 

MAJOR RUN CHARACTERISTICS Sccoild Dcltl^ H3t6 for Bottl LasgiT 3X1(1 USB Dopplsr 


O * « , IVT ^ 1 PGS-1040** 

Geopotentlal Model 

Lunar/Solar Gravitation __XES__ 
c —1 ^ 

Solar Radiation Parameter R 


Cd=2. 1 


Drag Parameters . 

Atmospheric Density Mode P«P»» Fil^l50 

„ , „ „ ^ State and Rho one 

Solve-For Parameters 


Editing Parameters 


3 Sigma 


USB-Doppler, Earth Tides 
Polar Motion, Marsh II 
Geodetics*** 


Arc 

Start 

Time 

Arc 

Length 

(hrs) 

No. 

of 

Sta- 

tions 

Observations 

Residual 

Statistics 

Maximum COMPARE 
Position Differences 
(m) 

Solve-For Parameters 
and Other Information 

Run 

ID 

Range 

Range-Rate 

No. 

Avail- 

able 

No. 

Used 

No. 

Avail- 

able 

No. 

Used 

Wtd. 

RMS 

standard 

Deviations 

Radial 

Cross- 

Track 

Along- 

Track 

Total 

RHO 

ONE 


PASSE 

5 

Range 

(m) 

Range- 
Rate 
(cm/ sec) 

780919 

30 

7* 


■ 



.83 


1.68 

-.65 


20* 



0.94 

11.76 

4.66 

12.28 

780920 

30 

9 



371 

325 

.99 


1.98 

-.67 


17 



1.01 

11.59 

2.15 

11.67 

780921 

30 

9 



366 

310 



1.93 

HH 


— 

20 



BRU 




780922 

30 

10 



513 

426 

.82 


1.64 

-.22 


25 



1.54 

3.21 

7.62 

7.70 

780923 

30 

9 


■ 

444 

392 

.82 


1.65 

-.21 


21 
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TABLE 5 


SATELLITE AND TIME PERIOD SEASAT - September 1978 

MAJOR RUN CHARACTERISTICS ^PPJ^oxlmately 30 Second Data Rate for Both Laser and USB Doppler 
mm !'GS-I040«« DT.B 3 sigma 

Lunar/Solar Gravitation Atmospheric Density Model H«P»> F#150 ntw Laser Range. Earth 

Solar Radiation Parameter Solve- For Parameters^_?55^® Tides, Polar Motion, Marsh II 

Geodetics*** 

















































































TABLE 6 


SATELLITE AND TIME PERIOD. 
MAJOR RUN CHARACTERISTICS 


SEAS AT - September 1978 

Approximately 30 Second Data Rate for Both Laser and USB Doppler 


Geopoteatial Model 


Lvinar/Solar Gravitation. 


PGS-1040** 


YES 


Solar Radiation Parameter. 


Cr=1.5 


Drag Parameters ^D ^ 


Atmospheric Density Modeftl 

State and Rho one 
Solve-For Parameters 


Editing Parameters ^ Sigma 

Laser Range and USB-Doppler , Earth 
Tides, Polar Motion, Marsh II 
Geodetics*** 


Arc 

Start 

Time 

Arc 

Length 

(hrs) 

No. 

of 

Sta- 

tions 

Observations 

Residual 

Statistics 

Maximum COMPARE 
Position Differences 
(m) 

Solve-For Parameters 
and Other Information 

Run 

ID 

Range 

Range-Rate 

No. 

Avail- 

able 

No. 

Used 

No. 

Avail- 

able 

No. ' 
Used 

Wtd. 

RMS 

standard 

Deviations 

Radial 

Cross- 

Track 

Along- 

Track 

Total 

RHO 

ONE 

■ 


■ 

Range 

(m) 

Range- 
Rate 
(cm/ sec) 

780919 

30 

6/7'* 

69 

66 

403 

345 

1.10 

1.52 

1.92 






0.67 

0.40 

2.22 

2.25 

780920 

30 

8/9 

79 

75 

371 

325 

1.15 

1.50 

2.06 

-0.71 


12/17 



0.93 

2.02 

3.25 

3.80 

780921 

30 

6/9 

89 

83 

366 

310 

1.16 

1.34 

2.08 

-0.4S 

IB 

14/20 



1.69 

0.91 

3.95 

4.00 

780922 

30 

8/lC 

79 

77 

513 

427 

0.99 

1.10 

1.84 

m 

■ 

m 



0.66 

2.72 

4. 18 

4.69 

780923 

30 

5/9 

64 

63 

444 

392 

0.96 

1. 11 

1.84 

HP 

B»IW» 

■■■■1 


10/21 
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TABLE 7 


SATELLITE AND TIME PEmnn SEASAT - September 1978 ^ 

ADDroximatelv Equal Laser and USB Doppler Observations Per Pass 
MAJOR RUN OHARAOTERrSTrCS ^PP^^^-^™^^^^-^^ ^Huct-L ^ 


Geopotential Model PGS-1040** 
Lunar/Solar Gravitation YES 

Cp=1.5 

Solar Radiation Parameter 


Cr,=2. 1 3 Sigma 

Drag Parameters n Editing Parameters 

,H.P. F#150 Laser Range and USB-Doppler, Earth 

Atmospheric Density Modeliil — L2 Qfher ^ . r;, c — t-i 

State and Rho one Tides, Polar Motion, Marsh ii 
Solve-For Parameters i ou*** 









































































































TABLE 8 


SATELLITE AND TIME PERTOD SEASAT - August 1978 

MAJOR RUN CHARACTERISTICS ^PP^^oximately 30 Second Data Rate for USB Doppler 


Geopotential Model 1040** 

Lunar/Solar Gravitation YES 

Cp=l 5 

Solar Radiation Parameter 


Drag Parameters 2 Editing Parameters ^ Sigma 

Atmospheric Density Model H.P., F#150 rn-TioT- USB— Doppler, Earth Tides 

Solve-For Parameters Polar MotiOIl, Marsh II 

(ieodetlcs*** “ 

















































































TABLE 9 


SATELLITE AND TIME PERIOD SEASAT - August 1 Q7« 

MAJOR RUN CHARACTERTSTTrs^PP^oximately 30 Second Data Rate for Laser 


Geopotential Model 1040** 
Lunar/Solar Gravitation 

Co=1.5 

Solar Radiation Parameter_U 


Drag Parameters — 5 Editing Parameters ^ Sigma 

Atmospheric Density Model Other ^Dge, Earth 

State and Rho one Tides, Polar Motion, Ms 



Arc Arc 

Start Length 

Time (hrs) 



780808 


780809 


780810 


30 6* 


30 6 152 130 


Observations 

Range 

Range 

-Rate 

No. 

Avail- 

able 

No. 

Used 

No. 

Avail- 

able 

No. 

Used 


Solve-For Parameters. 


Residual 

Statistics 


Maximum COMPARE 
Position Differences 
(m) 


Solve-For Parameters 
and Other Information 



Standard 

Deviations 


Range R^^ial ^ross- Mong- 


1 

780811 

30 

5 

142 

105 

780812 

30 

m 

105 

61 



2.05 

1.85 

2.03 

2.03 

1.56 

1.56 

2.42 

2.40 

1.99 

1.95 


Rate 
(cm/ sec) 



*Number ®f Stations 



















































































TABLE 10 


SATELLITE AND TIME PEmoD SEASAT - August 1978 

MA.TnppTTNPHARAPTERTSTrrs Approximately 30 Second Data for Both Lase r and USB Doppler 


Geopotential Model 


PGS-1040** 


Lunar/Solar Gravitation 

Solar Radiation Parameter. 




TV. . CtT=2. 1 

Drag Parameters ^ 

Atmospheric Density Mode^ 
Solve-For Parameters 


H.P., F#150 


3 Sigma 

Editing Parameters r 

Laser Range and USB-Doppler, Earth 


ides. Polar Motion, Marsh II 



Arc 

Arc 

Start 

Length 

Time 

(hrs) 



780808 


780809 


780810 


780811 


780812 


Observations 

Range 

Range 

-Rate 

No. 

Avail- 

able 

No. 

Used 

No. 

Avail- 

able 

No. 

Used 

^ 135 

83 

470 

400 

152 

122 

538 

429 

108 

102 

366 

317 

142 

103 

335 

276 

105 

60 

317 

269 


Residual 

Statistics 


Standard 

Deviations 


Range- 

Rate 

(cm/sec) 


Maximum COMPARE 
Position Differences 
(m) 


Solve-For Parameters 
and Other Information 


Radial 

Cross- 

Track 

1.53 

0.91 

0.66 

1.70 

1.36 

4.67 

2.30 

2.91 




































































































TABLE 11 


SATELLITE AND TIME PERTOD GEOS-III February and March 1976 

MAJOR RUN CHARACTERTSTTOS Approximately 30 Second Data Rate for USB Doppler 


Geopotential Model GEM lOB ** 
Lunar/Solar Gravitation 

Cr = 1.45 

Solar Radiation Parameter ^ 


Drag Parameters 

Atmospheric Density Model > E?!^75 

State Vector 
Solve-For Parameters 


Editing Parameters ^ Sigma 

ofw USB-Doppler, Earth Tides 
Polar Motion, Marsh II Geodetics*** 


! 

Arc 

Start 

Time 

Arc 

Length 

(hrs) 

760223 

30 

760224 

30 

760225 

30 

760226 

30 

760227 

30 

760228 

30 

760229 

30 

760301 

30 

760302 

30 


*Number of Stations and 
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TABLE 12 


SATELLITE AND TIME PERIOD GEOS-TII Fpbruary and Marr.h 1976 

MAJOR RUN OHARAOTERrsTTos Second Data Rate for USB Doppler 


Geopotential Model GEM lOB 
YES 

Lvmar/Solar Oravitatlnn 43 

Solar Radiation Parameter 


Drag Parameters 

H.P., F#75 

Atmospheric Density Mode^^^^^^ 

Solve-For Parameters 


Editing Parameters ^ Sigma 

TSB-Doppler, Earth Tides 

^5y a r Motion, Marsh II" ' Geodetic s *** 


Arc 

Start 

Time 

Epoch at 
15 hrs. 

Arc 

Length 

(hrs) 


Observations 

Residual 

Statistics 

Maximum COMPARE 
Position Differences 
(m) 

Solve-For Parameters 
and Other Information 

Run 

ED 

No. 

of 

Sta- 

tions 

Range 

Range-Rate 

No. 

Avail- 

able 

No. 

Used 

No. 

Avail- 

able 

No. 

Used 

Wtd. 

RMS 

Standard 

Deviations 

Radial 

Cross - 
Track 

Along- 

Track 

Total 

1 

■ 


1 

Range 

(m) 

Range- 
Rate 
(cm/ sec) 

760223 

30 

5* 


69, 

'159 i 

5/121 

.55 

0. 

8/1.2 



17* 



1.1 

10.8 

3.5 


760224 

30 

5 


94; 

187 ! 

8/138 

.51 

1. 

0/1.0 



21 



0.4 

5.2 

m 

5.4 

760225 

30 

5 


67, 

263 L 

6/193 

.60 

1. 

2/1.2 



25 



0.3 

2.3 

m 

2.5 

760226 

30 

4 


65, 

251 ' 

1/200 

.65 

1. 

1/1.4 



23 



HI 


M 


760227 

30 

4 


79, 

193 6 

7/163 

.63 

1. 

1/1.3 


■ 

21 



0.3 

m 

n 

8.2 

760228 

30 

5 


69) 

203 ! 

0/158 

.65 

1. 

7/1.2 



21 


■1 

1.5 

9.9 

BBli 


760229 

30 

5 


13) 

154 1 

1/110 

.61 

0. 

9/1.2 


■ 

■a 


■ 

0.3 

2.9 

2.2 

3.5 

760301 

30 

4 



.34 

114 

.57 


1.1 



19 



0.2 

5.2 

0.9 

5.3 

760302 

30 

5 


32) 

134 

8/105 

.72 

0. 

8/1.5 



15 

■ 




















HI 



m 

6.7 


*Number 

»f St£ 

tions 

and 

Passes 

1 for 

USB D 

jpplei 



■ 

■1 

■ 


■ 





**Computa 

lion 1 
1 IQ l 

ased 

Dn GE 
1 /f 

1 lOB 

- 1 /9< 

GM = 

ift 9R-; 

3986 

30.64 

km^/ se 

c^. 





■ 





^**Eilipso 

i / f 

d Pat 
= 1/: 

amete 

98.25 

*/ '■ 

rs fo 
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: Mars 

h 11 

Geode 

:ics : 

Rg=63) 
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TABLE 13 


SATELLITE AND TIME PERIOD GEOS-III Fphrnarv anH Marc-h 1Q7ft 

MAJOR RUN CHARACTERiSTirs Approximately 30 Second Data Rate for USB Doppler 


Geopotenttal Model lOB ** 
YES 

Lunar/Solar Gravitation 

Cg=T74o 

Solar Radiation Parameter 


Drag Parameters 5 ^ Editing Parameters ^ Sigma 

Atmospheric Density Mode l *^*> ^ rvftioT- Earth Tides 

s>tate vector ' Polar Motion, stun Geodetlcs*** 
Solve-For Parameters 


Arc 

Start 

Time 

Epoch at 
15 hrs. 

Arc 

Length 

(hrs) 

No. 

of 

Sta- 

tions 

Observations 


Residual 

Statistics 

Maximum COMPARE 
Position Differences 
(m) 

Solve- r or Parameters 
and Other Information 


Ra 

nge 

Range 

;-Rate 


No. 

Avail- 

able 

No. 

Used 

No. 

Avail- 

able 

No. 

Used 

wtd. 

RMS 

Standard 

Deviations 

Radial 

Cross- 

Track 

Along- 

Track 

Total 

1 

I 


1 

Run 

ID 

Range 

(m) 

Range- 
Rate 
(cm/ sec) 

760223 

30 

5* 


69 

'159 

'5/123 

1.88 

3 

1/4.0 


m 

17* 

wm 


5.3 

31.9 

an 

m 

760224 

30 

5 


94 

'187 

8/123 

1.08 

m 

0/2.2 



21 



10.4 

21.1 

sa 

9I 

■fiOSH 

760225 

30 

5 


67, 

'263 t 

6/186 

1.64 

■ 

2/3.1 



25 




15.0 

n 

15.3 

760226 

30 

4 


65> 

'251 f 

1/196 

1.90 

m 




23 

1 



in 


EH 

760227 

30 

4 


79y 

193 6 

2/157 

1.48 

3. 

2/3.6 

■ 


21 



3.2 

3.6 

HI 

n 

760228 

30 

5 


69> 

'203 / 

8/153 

1.59 

3. 

7/2.9 



21 



3.5 

14.5 

m 

99 

760229 

30 

5 


13, 

154 ] 

1/108 

1.50 

0. 

5/3.0 



14 



5.2 

29.9 

^9 


760301 

30 

4 



.34 

115 

1.21 


2.4 



12 



m 

9.2 

3.2 

9.5 

760302 

30 

5 


32, 

134 

8/98 

1.08 

1. 

8/2.1 



15 
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*Number ( 

StJ 

tions 

and 

Passes 

for 

USB D 

opplei 












**Computai 

:ion 1 
J39 1< 

ased 
m arid 

on GE 
1/f 
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= l/2( 
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8 . 9S7 
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30.64 

km^/se 
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= 1/2 
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TABLE 14 


Geopotential Model 


GEM lOB ** 


SATELLITE AND TIME p^Tiinn GEOS-lli March 1?7p 

MAJOR RUN rpAPArTPRTSTTUS Approximately 30 Second D ata Rate for U_SB Doppler 

** _ _ , Cn=3. 09 „ 3 Sigma 


Lunar/Solar Gravitation 

Solar Radiation Parameter 


Cr=1.45 


Cn=3.09 

Drag Parameters 

Atmospheric Density Model 1 

State Vector 

Solve-For Parameters. — — 


Editing Parameters r — 

U SB-Doppler, Earth Tides 
P^olar Motion, GEM 9 Geodetics*** 


Arc Arc 

Start Length 

Time (hrs) 

Epoch at 
15 hrs. 


760223 


760224 


760225 


760226 


760227 



Observations 

Range 

Range 

-Rate 

No. 

Avail- 

able 

No. 

Used 

No. 

Avail- 

able 

No. 

Used 



760228 


760229 


760301 


760302 
















































































































TABLE 15 


I 

OJ 

Ui 


SATELLITE AND TIME PERIOD GEOS-III February and March 197ft 


Geopotential Model 
Lunar/Solar Gravitation. 


MAJOR RUN CHARACTERiSTirs ^0 Second Data Rate for USB Doppler 
GEM lOB ** fT%=R DQ -3 0-* 

Drag Parameters D * t> ^ ^ Sxgma 


YES 


Solar Radiation Parameter 


Cr=1. 45 


HP F # 7 5 

Atmospheric Density Mode l * * * 

State Vector 
Solve-For Parameters 


Editing Parameters. 


USB-Doppler, Earth Tides 
^oTar Motion, Wgs Geodetlcs*** 


Arc 

Start 

Time 

Epoch at 
15 hrs. 


760223 


760224 


760225 


760226 


760227 


760228 


760229 


760301 


760302 


*Number 


***Ellipso 
H' 1/f 


Arc 

Length 

(hrs) 


30 


30 


30 


30 


30 


30 


30 


30 


30 


**Computa :ion based on GE| 
■ Raa=6378 1 139 hm and - 1/f 


id 


Pait'' 


No. 

of 

Sta- 

tions 


5* 


5 


pf Stations 


Observations 


Range 


No. 

Avail- 

able 


No. 

Used 


69 


94 


[187 58/13a 


and 


ametelrs fop: 



67 


'263 46/194^ .85 


65 


79 


69 


13 


32 


Passes 


;M lOB 


WGS 


Range-Rate 


No. 

Avail- 

able 


159 55/117^ 


251 31/199 


193 47/163 


h54 


134 


134 


for 


GM 4 


Geode 


No. 

Used 


Residual 

Statistics 


Wtd, 

RMS 


.43 


[203 51/153 




114 


8/103 


tics : 


.75 


.77 


,78 


76 


,84 


standard 

Deviations 


Range 

(m) 


8/0.9 


,55 


,68 


USB Dbpplef 


398600.64 


0 


km^/ 


se 


Rp=6B78.13$ 


Range- 

Rate 


2/1.5 


7/1.6 


5/1.5 


6/1.4 


9/1.2 


4/1.5 


1.1 


6/1.4 


c2, 


km 


Maximum COMPARE 
Position Differences 


Radial 

Cross- 

Along- 

Total 


Track 

Track 


IB 

20.8 

6.5 

mm 

0.4 

m 

1.8 

5.0 

IB 

8.0 

3.6 

8.2 

0.7 

9.2 

2.8 

9.5 

0.5 

6.0 

3.3 

6.7 


13.4 

n 


'IB 

6. 1 

5.6 

7.7 

0.6 

12.6 


n 




HI 




HI 




H 














Solve-For Parameters 
and Other Information 


PASSES 


17* 


21 


25 


23 


21 


21 


14 


12 


15 


Run 

ID 


ani 
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TECHNIQUES FOR INCREASING THE EFFICIENCY OF EARTH GPAVTTY 
CALCULATIONS FOR PRECISION ORBIT DETERMINATION * 

Richard L. Smith 
Anatoly S. Lyubomirsky 
Computer Sciences Corporation, 

ABSTRACT 

Two techniques for increasing the efficiency of Earth grav- 
ity calculations are analyzed. The first is a representa- 
tion using Chebyshev expansions in three-dimensional cells. 
Mathematical formulas are given for converting the standard 
spherical harmonic representation fe.g., GEMlOB 36 x 36) to 
the Chebyshev representation. The error in the truncated 
Chebyshev representation was measured as a function of cell 
size and degree of truncation. For example, with a sixth 
degree Chebyshev expansion, the maximum gravity error is 
about 10 ^^g for a 36 x 36 parent representation in a cell 
extending 5 degrees in both latitude and longitude and hav- 
ing a thickness of 600 kilometers. Computer storage re- 
quirements and relative CPU time requirements are presented. 
The Chebyshev gravity representation can provide a signif- 
icant reduction in CPU time in precision orbit calculations, 
but at the cost of a large amount of direct-access storage 
space, which is required for a global model. 

The second technique employs a temporary file for storing 
the components of the nonspherical gravity force. In 


*This work was supported by the Operations Analysis Section, 
Operational Orbit Support Branch, Goddard Space Flight 
center. National Aeronautics and Space Administration, under 
Contract NAS 5-24300. 
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differential correction orbit solutions it is often unneces- 
sary to repeat computations for most of the gravity terms 
during subsequent iterations for which the satellite's posi- 
tion changes only slightly. By saving a direct-access file 
of gravitational forces and partial derivatives it is pos- 
sible to reduce CPU time without significantly affecting 
orbit accuracy. The gravity file is updated whenever the 
position tolerance is exceeded. The Goddard Trajectory De- 
termination System was temporarily modified to test this 
technique, and the results of the test are presented. 

1. INTRODUCTION 

As the orbit determination accuracy for Earth-orbiting 
spacecraft is improved through the use of increasingly more 
accurate Earth gravity models, the computer time require- 
ments increase rapidly. Using the customary global spheri- 
cal harmonic expansion, the amount of computation time 
increases approximately as the square of the maximum degree 
and order of the expansion. For currently available gravity 
models, for example, the Goddard Earth Model lOB (GEMlOB) , 
most of the computation for an orbit solution is devoted to 
evaluations of the gravity force. Clearly, less time- 
consuming methods of gravity evaluation are required, par- 
ticularly if precise gravity models are needed for future 
operational orbit determination. The need for faster meth- 
ods is enhanced by the fact that the utilization of more 
precise gravity models requires the use of correspondingly 
smaller step sizes for numerical integration of the space- 
craft equations of motion. 

Table 1 shows the amounts of computer time (GSFC IBM 
S-360/75) currently required for orbit solutions calculated 
using the Goddard Trajectory Determination System (GTDS). 

In order to isolate the dependence of the computer time on 
the specified value of the maximum degree and order in the 
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Table 1. GTDS Computer Time Usage for Various Sizes of 
the Spherical Harmonic Gravity Expansion 

SPACECRAFT: SEASAT-1 

NUMERICAL INTEGRATOR: COWELL FIXED STEP, 12TH ORDER 
FORCE MODEL: 

• GRAVITY: SOLAR, LUNAR, GEM9 

• DRAG, WITH HARRIS-PRIESTER ATMOSPHERE 

• SOLAR RADIATION FORCE 

• MEAN OF 1950,0 SYSTEM FOR INTEGRATION 

EPOCH: 18*^ ON JULY 10, 1978 ARC LENGTH: 30 HOURS 
EPOCH - ARC LENGTH: 18*^ ON JULY 10, 1978 - 30 HOURS 
OBSERVATIONS: 391 DOPPLER USB, 100 LASER RANGE 



IBM S-360/75 COMPUTER TIME USAGE (MIN) 

GRAVITY MODEL 

90-SECOND STEP SIZE 

45-SECOND STEP SIZE 


CPU 

I/O 

CPU 

I/O 



EPHEM PROGRAM 


4x4 

0,888 

0.241 

1.544 

0.239 

8x8 

1.007 

0.241 

1.613 

0.239 

21 X 21 

1,280 

0.252 

2.306 

0.249 

36 X 36 (GEM10B) 

3.210 

0.329 

5,058 

0.330 



DC program' 


4x4 

7.448 

1.804 

11.015 

1.725 

8x8 

8.322 

1.805 

12.051 

1.727 

21 X 21 

10,419 

1.817 

15.482 

1 .739 

36 X 36 (GEM10B) 

20.577 

1.938 

35.952 

1.855 


'six ITERATIONS AND CONVERGENCE 
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spherical harmonic expansion, all other input parameters for 
these solutions were identical. Computer times for both GTDS 
Ephemeris Generation (EPHEM) and GTDS Differential Correction 
(DC) Program runs are shown in this table. 

TWO methods for efficiency improvement are examined in this 
paper. Section 2 outlines a gravity representation using 
Chebyshev polynomials rather than spherical harmonics. Sec- 
tion 3 considers a procedure for making use of previously 
computed values of the gravity force during the later itera- 
tions of differential correction orbit solutions. This 
procedure, unlike the Chebyshev representation, is not gen- 
erally applicable to orbit prediction. Section 4 assesses 
the merit of these two methods and indicates directions for 
future work. 

2. REPRESENTATION OF THE EARTH'S GRAVITY FIELD 
USING CHEBYSHEV POLYNOMIALS 

2.1 OUTLINE OF THE METHOD 

in order to accurately represent the Earth's gravity using 
Chebyshev polynomials, the region of interest is partitioned 
into cells, and for each cell the gravity force components 
are expressed as a series of Chebyshev polynomials. The 
numerical values of the expansion coefficients for a given 
cell are, in general, different from those of any other 
cell. With a suitable selection of the cell dimensions, the 
convergence of the Chebyshev series is sufficiently fast 
that the computational effort for its evaluation is signifi- 
cantly less than the effort required to evaluate the stand- 
ard spherical harmonic expansion. In exchange for the 
reduction in computational effort, however, the Chebyshev 
representation requires a large data set containing the ex- 
pansion coefficients for all of the cells. 
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The evaluation of the gravity, force in GTDS is accomplished 
with the following standard spherical harmonic expansion: 


n 


max 




n 


n 


■g} (n + 1 ) I ^ fsin (i>) 
n=0 m=0 


( 1 ) 


cos mX + sin mX 
n n 


max 


= 9 


E (^)" E 


(sin ({)) - m tan (j) (sin (j>) 


n=0 


m=0 


M2) 


cos mX + sin mX 
n n 



g 

cos (j> 


n 

max 

E 

n=0 



n 


E 

m=0 


m P 


( sin 


4 >) 


cos mX 


sin mxj 


( 3 ) 


where r = radial distance in Earth radii (a) 

(|) = geocentric latitude 
X = geocentric longitude 

p™ = Legendre function of degree n and order m 

rimax “ maximum degree of the spherical harmonic expan- 
sion for the Earth's gravity field 

g = GM/(ar)2, where G is the universal constant of 
gravitation, M is the Earth's mass, a is the 
Earth's radius, and r is defined above 


C'^ = nonnormal ized spherical harmonic expansion co- 

efficients for the geopotential field model con- 
sidered 
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The Chebyshev expansions used in this paper, also yield the 
radial, latitudinal, and longitudinal gravity components, F^, 
F^, F^. The Chebyshev expansions are applied only to that part 
of the gravity force described by spherical harmonic terms 
of degree greater than 4. Terms of degree less than or equal 
to 4 are still evaluated using spherical harmonics. 

In each cell, independent position variables, x, y, and z, 
are designated. These variables are related to r, (j), and X 
by means of the following equations: 


|=^+Ax (lxl<l) (4) 

o 

sin ()) = sin 4)^ + Cy ( I (}) I <_ 45°, lyl £ 1) (5) 

cos (j) = cos + Cy ( I (}> I > 45°, lyl £ 1) (6) 

cos X = cos X^ + Dz ( IX - 90° 1 <_ 45°, Iz I <_ 1) (7) 


The cell origin is (r , d) , X ) and the phys ical size of a 
cell is controlled by the three parameters A, C, and D. The 
position variables x, y, and z describe displacements, rela- 
tive to the cell origin, in the radial, latitudinal, and 
longitudinal directions, respectively. The locus of points 
such that X = +1 or X = -1 describes spherical surfaces 
bounding the top and bottom of a cell. The locus of points 
such that y = +1 defines cones of constant latitude bounding 
the north and south sides, and the locus of points such that 
X = +1 describes longitudinal planes bounding the cell on 
the east and west sides. This choice of independent vari- 
ables leads to cell crowding near the poles, but allows a 
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fast and efficient procedure for calculation of the Chebyshev 
expansion coefficients. 

AS indicated by Equations (5) and (6) , the latitude-like 
variable, y, is defined differently for the polar and equa- 
torial regions. This difference is necessary to avoid slow 
convergence of the Chebyshev expansions close to the poles 
and close to the equator. This slow convergence problem also 
exists for X = 0 or X = tt using the definition given for z 
by Equation (7). However, it is only necessary to apply a 
longitude shift when the problem occurs (by suitably adjust- 
ing the cpj's and S{^'s) and thus avoid a double definition. 

The expansion of each factor of a typical spherical harmonic 


_1 

n+1 ^n 
r 


( sin 


4 >) 


cos 

sin 


mX 


into a series of Chebyshev polynomials follows the equations 
(for each cell) 




k=0 


2 



Z 


( 1 ) 

mk 


Tk(z) 


( 8 ) 


(9) 


cos mX 


7T 


( 10 ) 



sin mX 


-z 

k=0 


'2 - 5, 


<0 


Z ^ ^ ^ T ( Z ) 
^mk ^k'^^ 


( 11 ) 


The Chebyshev polynomials, , are functions of x, y, or z and 
satisfy the recurrence relation 


Ti_|_i(x) = 2x T^(x) - T^_^(x) 


( 12 ) 


where the subscript indicates the degree of the polynomial. 

In several cases, the Chebyshev expansions indicated by Equa- 
tions (8) through (11) are finite, not infinite, as a result 
of the definitions of x, y, and z. The X's, Y's and Z's are 
the Chebyshev expansion coefficients and their values depend 
on the cell parameters r^, 4>^, X^, A, C, and D, in addition 
to the order and degree of the spherical harmonic. 

The X's, Y's, and Z's are combined in the following way, 
according to Equations (1) through (3), to form the three 
subscripted Chebyshev expansion coefficients, e.g., 
used for the calculation of the force components: 

ma x n 

= « Z ^nl Z (o 

n=4 m=0 


+ S^Z 

mk n mk 


(13) 



‘‘‘ma x h 

0 Z Z 

n=4 m=0 


m+1 

nj 


n mk 


+ 


S 


( 2 )\ 

n'^mk / 


(14) 


n 


max 


n 


c ; 


(3) 


ijk 


= Q 


Z >..L 


mY 


m 

nj 


n=4 


m=0 


n mk 


( 2 ) 

^n^mk 


(15) 
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max 


C 


(4) 


ijk 


== Q 


L E 


mY 


m 

nj 


n=4 


m=0 


s”z 

n mk 


C 


(2)\ 
n'^mk j 


(16) 


(2 - 6 


Q = 


Oi^ 


(2 - 6 q .) (2 


'^Ok^ 


(17) 


The three gravity force components are then calculated in 
the following way: 


IJK. 

f'r = -9 E E E ‘^Ijk '18) 

i=0 j=0 k=0 
IJK 

% ° 9 5] E E Hjk - Ti(x)T.(y)T|^(z) (19) 

i=0 j=0 k=0 ^ 

IJK 

I'l = ssfmF E E E ^ijk '8“> 

i=0 j=0 k=0 

These three equations represent the calculation of gravity 
as it might be performed in an orbit determination program, 
using precalculated coefficients. 

The formulation used in this paper required four types of 
three-subscripted Chebyshev expansion coefficients. With 
additional work, it should be possible to also expand the 
function 


tan (j) (sin (J)) 
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in a Chebyshev series, leading to a formulation using only 
three types of coefficients. This additional complication 
was omitted for the present for simplicity. 

AS indicated by Equations (8) through (16) the three- 
subscripted coefficients depend on the gravity model coeffi- 
cients, and s|^, the cell location, and the cell dimensions. 
The combined set of three-subscripted coefficients for all 
cells constitutes a Chebyshev representation for the given 
gravity model. 

The calculation of the Chebyshev coefficients for the spher- 
ical harmonic factors, that is, the calculation of the X's, 
Y's, and Z's, can be easily accomplished using recurrence 
relations. These recurrence relations are as follows: 

Recurrence relations for the radial Chebyshev coefficients: 


^n + l,i 2 ^^n,i+l ^n,i-l^ ''' r^ ^n,i 


(n > 0 , all i) (21) 


X 


1 


n+1,0 n 4- 1 


='n,0 


- - y Vl,0 


( 22 ) 


(n > 0) 


Recurrence relations for the longitudinal Chebyshev coeffi- 
cients: 


2mil,k = ^ ^ ^mHtl) 'O 

- Z , (all m, all k ) 

m “ i. / K 


(23) 
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„ ( 2 ) ^ ^ / ( 2 ) , 7 ^ 2 ) \ ^ ^ ( 2 ) 
ni+i,i< ^ m,k-l m,k + l/ 0 m,k 


( 24 ) 


Recurrence relations for the latitudinal Chebyshev coeffi- 
cients ( 1 (|) i £ 45°) : , 




'n+1 , j 


2n + 1 C m m 

n - m 4- 1 2 \ n,j-l n,j+l 


TT ^0 ^n,j 


2n + 1 
n - m 

n + m 
n - m + T ■"n-l,j 


m 


yin 


(all j, n >m > 0) 


(25) 


■n, j 


(2n - 1) (2n - 3) 

,n-2 


4 


C sin (})q ^n-2,j+l 


n-2 n-2 

^n-2,j-2 ^n-2,j+2 


,n- 2 


(26) 


+ ( 1 - sin 4). 


Cl\ yn-2 
2 / ^n-2,j 


(all j, n > 2) 
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Recurrence relations for the lati'tudinal Chebyshev coeffi- 
cients ( I (j) I > 4 5° ) : 


) 

n+2,i (n + 1 - m) (n + 2 - m) | 


(2n + 3) 


+ 


- (2n + 1) (C cos 4)^) 


(n + m)(n + iT\-l) 
(2n - 1) 

m m 

^n,i+l ^n,i-l 




n- 2 , i 


m 




- '2" - n I-Jl^n.l-^2 + '^n,l-2 

(n + 1 - m) (n + 1 + m) 

(2n + 3) 


(27) 


+ (2n + 1) sin'^ ~ (2n + 1) y 


(n + m) (n - m) 


(2n - 1) 


n, 1 


( all i , n > m > 0 ) 


Cli = (2" - n 


+ Y 


n 


n , i -1 


. „n , C /,,n 

*0 ^n,l * 2 Pn,i+1 


(all i , n > 0 ) 


(28) 


The derivation of these recurrence relations is omitted 
here; some detail is given in Reference 1. It should be 
noted that, although the same symbol is used in each case, 
the Y's of Equations (25) and (26) are defined differently 
than the Y's of Equations (27) and (28). There should be no 
confusion since Equations (25) and (26) are intended only 
for the equatorial region, while Equations (27) and (28) 
apply to the polar regions. 


2 . 2 ERROR MEASUREMENTS FOR THE CHEBYSHEV REPRESENTATION 

This section addresses the question of how closely a 
Chebyshev gravity representation matches the gravity field 
defined by the parent spherical harmonic representation. In 
order to study the Chebyshev expansion error, a computer 
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program was written to numerically evaluate the error for 
any selected cell. The program first constructs the 
Chebyshev expansion coefficients for the given spherical 
harmonic expansion, using the recurrence relations given in 
Section 2.1. These Chebyshev expansion coefficients are 

functions of the C^'s and S^'s; the cell parameters r^, 
snd and A, C, and D. Then, for a selected maximum 

degree, the three gravity force components, F , F, and F, 

r (p A 

generated by the Chebyshev expansions (Equations (18) 
through (20)) are numerically compared with the corresponding 
force components calculated from the spherical harmonic ex- 
pansion (Equations (1) through (3)), using a minimum degree 
of 4. This comparison is made at many points uniformly . 
distributed throughout the given cell, and the maximum dif- 
ference between the two representations provides a measure 
of the Chebyshev expansion error. All of the error measure- 
ments in this paper apply to Chebyshev representations based 
upon the GEMlOB 36 x 36 gravity model. 

Figures 1 and 2 show the numerically computed error as a 
function of the cell size parameter A. For simplicity, the 
latitude size parameter C, and the longitude size param- 
eter D, remained equal to A as A was varied. Figures 1 and 
2 show the error for cells at reference heights of 967 kilo- 
meters and 255 kilometers, respectively. On each figure, a 
reference error level at 10''^^g is indicated. Order of 
magnitude estimates place the resultant orbit error at less 
than 0.1 meters for a 5-day orbit propogation subject to a 
high-frequency gravity error having this amplitude. The 
maximum degrees for each of the Chebyshev components were 
equal to one another and are indicated for each group of 
curves in the figure. For example, in Figure 1, the upper 
group of curves represents the error in the three-force com- 
ponents as a function of A for a 3 x 3 x 3 Chebyshev expan- 
sion. 
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CHEBYSHEV EXPANSION ERROR (IN UNITS OF 9) 



Figure 1. Numerical Measurement of Chebyshev Gravity 
Representation Error as a Function of Cell 
Size and Expansion Degrees (Height of Cell 
Center = 967 Kilometers) 
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CHEBYSHEV EXPANSION ERROR (IN UNITS OF g) 



0.001 0.002 0.01 0.02 0.1 


CELL SIZE PARAMETER A (A = C = D) 


Figure 2. Numerical Measurement of Chebyshev Gravity 
Representation Error as a Function of Cell 
Size and Expansion Degrees (Height of Cell 
Center = 255 Kilometers) 
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Each of the error curves in Figures 1 and 2 has a range, for 

intermediate values of A, where the curve is nearly a 

straight line. In this range, the slope of this straight 

line, on a log-log scale, is one greater than the maximum 

degree of the Chebyshev expansion; i.e., the error varies as 

the cell size to the power, where is the 

maximum Chebyshev degree. (This rule does not seem to be 

accurate for the larger values of For larger 

values of A, the curves bend away from the straight line. 

For very small values of A, a numerical noise level is 

- 1 8 

reached and the error reaches a lower limit--about 10 g 
for Figure 1 and 3 x 10^"^g for Figure 2. 

Figures 3 and 4 show the numerical error as a function of 
latitude for a 5° x 5° cell, using a 6 x 6 x 6 poly- 
nomial degree expansion. The cell thickness was chosen to 
be small, at a value of 12.8 kilometers, to eliminate the 
effects of radial variation on the error. The results in 
Figure 3 were obtained using the equatorial zone formulation 
(Equations (5), (25), and (26)) and those in Figure 4 were 

obtained using the polar zone formulation (Equations (6) , 

(27) , and (28)). The former diverges near the poles and the 
latter diverges near the equator, so that a global Chebyshev 
gravity model must be based upon a combination of these two 
formulations, in Figures 3 and 4, the maximum error in each 
cell is plotted at the cell center, so that cells centered 
at 2.5 degrees latitude extend to the equator and cells cen- 
tered at 87.5 degrees extend to within 0.001 degrees of the 
pole. 

The slight rise in error near the pole in Figure 4 occurs at 
error sampling points that are 0.75 degrees from the pole. 
This slight rise is presumably due to factors of cos 4) and 
an associated loss of precision in the calculation of F^ 
and F^ (Equations (2) and (3)). 
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CHEBYSHEV EXPANSION ERROR (IN UNITS OF 
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CHEBYSHEV EXPANSION ERROR (IN UNITS OF g) 




GEM10B 36 X 36 

AX -5 DEGREES 

A0 = 5 DEGREES 

Ar -12.8 KILOMETERS 

MAXIMUM CHEBYSHEV DEGREE :6 

r -7347 KILOMETERS 
o 

X -32.6 DEGREES 



LATITUDE (DEGREES) 


ERROR IN f. 


ERROR IN F 
ERROR IN F, 


10 ' 



•P 


Figure 4. Numerical Measurement of Chebyshev Gravity Representation 

Error As A Function of Latitude (Polar Zone Expansion Used) 
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Outside of the latitude regions in which divergence of the 
Chebyshev expansions is approached, it is clear from Fig- 
ures 3 and 4 that a uniform level of error is obtained using 
cells of constant latitudinal and longitudinal dimensions. 
The solid angle of these cells is much smaller near the 
poles than near the equator; leading to an unpleasant crowd- 
ing of cells near the poles in a global Chebyshev model. 

2*3 ESTIMATED CHARACTERISTICS OF A GLOBAL CHE BYSHEV GRAVITY 
REPRESENTATION ^ 

The use of the Chebyshev representation for precise satel- 
lite orbit determination requires a large, direct-access 
data set that contains the three-subscripted Chebyshev coef- 
ficients for a distribution of cells covering the entire 
spatial region of interest. The orbit determination program 
would retain in main memory the coefficients for a small 
number of cells and would update this working storage as 
necessary, drawing from the large, direct-access data set. 

In this section the general characteristics of a sample 
global Chebyshev representation are estimated. 

Table 2 provides data for estimating the speed of the 
Chebyshev representation, relative to the spherical harmonic 
representation. For each representation, the table shows 
the number of machine multiplication or division operations 
required to evaluate the three force components at a single 
spatial point. The numbers given assume efficient coding. 
The maximum degree used in the Chebyshev representation, 

^max' assumed to be chosen to be the same for all three 
indices in the expansions. Comparing the 36 x 36 spherical 
harmonic representation with the 6x6x6 Chebyshev repre- 
sentation, the latter requires about 75 percent less time 
for force evaluation (1,736 operations versus 6,933 opera- 
tions) . . 
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Table 2. Number of Computer Multiplication or Division 
Operations Needed for Gravity Force Evalua- 
tion in the Chebyshev and Spherical Harmonic 
Gravity Force Representations 


CHEBYSHEV REPRESENTATION 

MAXIMUM DEGREE 

NUMBER (Ni) 

OF MULTIPLICATIONS 
OR DIVISIONS" 

3 

332 

4 

640 

5 

1,098 

6 

1,736 

8 

3,669 

10 

6,685 


*Ni =5(Kmax + 1)3 + 3Kmax 


SPHERICAL HARMONIC REPRESENTATION 

MAXIMUM DEGREE 
I'^max* 

NUMBER (N2) 

OF MULTIPLICATIONS 
OR DIVISIONS"" 

, 4 

1 16 

8 

409 

16 

1,473 

21 

2,463 

30 

4,875 

36 

6,933 

48 

12,129 


*"*M2 - 15 
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since the number of operations in the Chebyshev representa- 
tion increases as the third power of K , while the num- 

max 

ber of operations in the spherical harmonic representation 

increases as only the square of the maximum degree, it is 

desirable to choose as small a value as possible for K 

■ max 

in order to achieve a computation time advantage. In order 
to simultaneously meet accuracy requirements, it is then 
necessary to properly adjust the cell dimensions. 

The characteristics of the Chebyshev model presented in Fig- 
ure 5 were based upon Table 2 and the results of Sec- 
tion 2.2. This sample model covers the range of many NASA 
low-altitude spacecraft; an additional layer could be added 
to extend the model to higher altitudes. The estimate of 
the total number of three-subscripted Chebyshev coefficients 
assumes that only three types were necessary. Although the 
formulation presented in Section 2.1 employed four types of 
these coefficients, it is expected that there would be no 
difficulty in modifying the formulation to require only 
three types. 

From Figure 5, it is clear that the computation time advan- 
tage of the Chebyshev representation is accompanied by the 
need for a large, but not unreasonable, amount of direct- 
access storage. 

3. FILE RETRIEVAL FOR GRAVITY FORCE EVALUATION 
3 . 1 FILE RETRIEVAL METHOD 

In standard GTDS Differential Correction orbit solutions, 
the full force model is reevaluated during every iteration. 
Except for the first and second iterations, corrections to 
the orbital position are generally so small that the change 
in position has a negligible effect on the numerical values 
of most of the spherical harmonic terms in the gravity model. 
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• ACCURACY: 10“^°g FOR GEM10B 36 X 36 

• MAXIMUM DEGREE OF EXPANSION: 6 X 6 X 6 

• NUMBER OF CHEBYSHEV 

COEFFICIENTS FOR EACH CELL: 3 x (7 x 7 x 7) = 1029 

• CELL SIZE; Ah = 607 Kl LOMETERS (A = 0.04) 

A0 = 5 DEGREES 
A\ = 5 DEGREES 

• CELL DISTRIBUTION; SINGLE LAYER (r^ = 6954 KILOMETERS) 

hMiN = 284 KILOMETERS 

• NUMBER OF CELLS: 36 X 72 = 2592 

• NUMBER OF CHEBYSHEV 

COEFFICIENTS IN STORAGE; 2592 x 1029 = 2.7 MILLION 

• CPU TIME FOR GRAVITY EVALUATION 
(RELATIVE TO SPHERICAL HARMONICS): 0.25 


Figure 5. Characteristics of a Sample Chebyshev 
Gravity Model 
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Rough estimates have indicated that, for a 1-day orbit, a 
10-meter error in the argument of the portion of the gravity 
force that does not include the monopole and quadrupole 
terms leads to orbital position errors that are well below 
0.01 meter. These estimates suggest that considerable compu- 
tation time could be saved, particularly for a 36 x 36 grav- 
ity model, if a file of gravity values was saved for use 
during the later iterations. 

The method of gravity evaluation tested is shown in Figure 6. 
This figure is a flowchart representing the GTDS subroutine 
that evaluates the gravity force, F(N x N) , for a given in- 
put position. A test is first made to determine whether a 
gravity file value exists for the given integration point. 
(This method is valid only for fixed-step numerical integra- 
tion.) If the file value exists, then the position associ- 
ated with the file is compared with the input position. If 
the difference is less than a prescribed tolerance, e, 
then the file value is accepted. The file value describes 
that part of the gravity force represented by spherical har- 
monic terms of degree greater than four. This value is ad- 
ded to the 4x4 force calculated for the input position, 

F(4 X 4) , to produce the total gravity force F(N x N) . 

If the file gravity value does not exist, or if the position 
deviation 1^1 is greater than the specified tolerance, s, 
then the file is not used. Instead F(N x N) , F(4 x 4), and 
F(FILE) are calculated, F(FILE) is stored for later use, and 
F(N X N) is returned by the subroutine. The resultant orbit 
precision of this method is controlled by the specified 
value of e. 

Not shown in Figure 6 is the treatment for partial deriva- 
tives of the gravity force with respect to position. These 
are stored, retrieved, and calculated in a manner parallel 
to that of the force components themselves. 
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Figure 6. Method for Gravity Force Evaluation Using 
File Retrieval 
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3 . 2 FILE RETRIEVAL RESULTS 

in order to test the file retrieval method, two GTDS differ- 
ential correction orbit solutions, 12 hours in length, were 
calculated using a 36 x 36 Earth gravity model and using 
Unified S-Band and laser tracking data. One solution was 
calculated in the standard way, and the other used the file 
retrieval method. For the latter solution, the position 
tolerance, e, was specified to be 500 meters. Each solu- 
tion required four iterations to converge, and each differ- 
ential correction solution was followed by 12-hour ephemeris 
generation, using the converged orbital elements. The a 
priori elements for the two solutions were identical, dif- 
fering from the converged elements by about 80 meters. 

A direct comparison between the ephemerides of the two solu- 
tions is shown in Figure 7. The position difference between 
the two solutions is plotted over the solution time inter- 
val. Examination of the intermediate results showed that 
for the first hour, the gravity file was built, but never 
subsequently updated since the 500-meter tolerance was never 
exceeded. On the other hand, for the following 11 hours, 
the gravity file was built during the first iteration, and 
since the 500-meter tolerance was exceeded during the second 
iteration (because the first-iteration orbit error progres- 
sively worsened with time, and this first-iteration orbit 
was the basis for the first-iteration file) the file was 
automatically updated, using positions generally accurate to 
5 meters. The last two iterations were calculated with no 
further updates to the file. This file update history ex- 
plains the sharp drop in orbit error over the first half 
hour in Figure 12--from 42 millimeters to the 5-millimeter 
level. 

It is clear from this file update history that the file re- 
trieval method reduces the number of standard gravity force 
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POSITION ERROR (MILLIMETERS) 


ORIGINAL FILE 



Figure 7. Orbit Error Resulting From Use of Gravity File With Position 
Tolerance Specified at 500 Meters 
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evaluations by more than a factor of two without substantial 
orbit precision loss. The CPU times for the two solutions 
were 1.23 minutes and 0.69 minutes (IBM S-360/95) for the 
standard and file retrieval solutions, respectively. These 
CPU times do not accurately show the full potential computa- 
tion time reduction of the file retrieval method because, 
for simplicity, these test calculations did not incorporate 
file usage into the numerical integration starting algo- 
r ithms . 


4. CONCLUSIONS 

The results presented in this paper show that the Chebyshev 
representation should provide substantial computation time 
savings for orbit determination using precise Earth gravity 
models, although its disadvantage is the requirement for a 
large file of pre-calculated Chebyshev coefficients. Tests 
of this representation in actual orbit calculations need yet 
to be performed. 

Two areas for possible improvement for the Chebyshev repre- 
sentation are evident. First, truncation of terms in the 
three-dimensional expansion should be explored. Rather than 
summing over terms such that i, j, and k range from 0 to 
K , it may be possible to sum over terms such that 

ITla, X 

i + j + k ranges from 0 to K . This tvpe of summation 

max * 

reduction could save a factor of approximately three in both 
execution time and in direct-access storage. The second 
improvement would be to extend the formulation so that 
Cartesian components of the gravity force are directly cal- 
culated, rather than spherical components. This would re- 
quire the derivation of additional recurrence relations for 
evaluation of the Chebyshev coefficients. 

The file retrieval method for gravity evaluation has been 
shown to be an effective method for reducing computation . 
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time without sacrificing orbit accuracy. 
Chebyshev representation, it could almost 
tion time problems in orbit determination 
available, precise gravity models. 


Combined with the 
eliminate computa 
using currently 
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AN ANALYSIS OF SIMULTANEOUS SATELLITE VISIBILITY 
TIME SPANS FOR TWO EARTH OBSERVATION STATIONS 

F. K. Chan 
Phoenix Corporation 

1700 Old Meadow Road, McLean, Va. 22102 
ABSTRACT 

Analysis was performed to estimate the statistical 
visibility time spans of earth orbiting satellites as seen 
simultaneously by a groiand station and a ship. The analysis 
covers topics such as time average population, average population 
times and also the percentage visibility times for a given 
number of satellites. These results are useful for specific 
communications satellite applications. Numerical results are 
obtained for various configurations of ground station and ship. 
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SECTION 1 


INTRODUCTION 


This report is concerned with the analysis of the number and also 
the time of satellites mutually observed by both a ground station and a 

ship. Unlike the relatively simple case of a single observation sta- 
tion for which the region of observation is the volume bounded by a 
cone, the present more' complicated case has a region of observation de- 
termined by the intersection of two cones. This region has a volume 
determined only by the separation distance between the ground station 
and the ship; but it also has a directional property determined by the 
relative position of the ship with respect to the ground station. 

Because the analysis becomes extremely complex, it is necessary to make 
certain simplifying assumptions. 

The first assumption is that the satellites presently orbiting the 
earth may be broadly classified into a few categories. This simplifica- 
tion is supported by the fact that^^^ since 1977 approximately 635 sat- 
ellites have been launched and these may be characterized as in Table 1.1 




Table 1.1 



Class 

J^verage 

Period 

Average 

Inclination 

Average 

Altitude 

Number 

I 

100 min. 

O 

O 

00 

800 km 

440 

II 

12 hr. 

o 

O 

20,000 km 

106 

III 

24 hr. 

0° 

36,000 km 

57 

IV 

Others 



32 

Thus , 

instead of having 

to deal with the volume 

of the region 

of observe 


tion, the analysis deals with the areas at the various altitudes. In 
this analysis, only Class I and II satellites are considered. Class III 
satellites are considered separately because they are geosynchronous. 
Class IV satellites are irregular and will not be considered at all. 


(1) NASA, Satellite Situation Report, Volume 21, Number 1, 
February 28, 1981. 
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The second assumption is that within each of the two categories con 
sidered, the satellites have circular 'orbits which are uniformly distrib 
uted in terms of equatorial crossing and, moreover, the satellites are 
also uniformly distributed along the orbital arcs. 

Section 2 deals with the derivation of the number density of sate- 
llites in this statistical distribution. Section 3 deals with the deter 
mination of the common region of observation of both the ground station 
and the ship. Section 4 is concerned with the computation of the time 
average population of satellites within the mutual region of observation 
Section 5 briefly discusses the computation of the average population 
times of these satellites in the same region. Section 6 summarizes the 
results of this sfudy for Class I and II satellites. 

Readers who are strictly interested in the numerical results may 
go directly to Section 6 and omit the intervening sections which deal 
with the mathematical analysis. 
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SECTION 2 


STATISTICAL DESCRIPTION OF ORBITING SATELLITES 


2 . 1 Distribution Function 

Consider a statistical description of a system of N satellites as 
previously described in which the circular orbits are uniformly distributed 
in terms of equatorial crossing and the satellites are uniformly dis- 
tributed along the orbital arc. Consider Figure 2 . 1 which illustrates 
a given orbit with inclination i. Let 0 be the latitude, (j) be the right 
ascension measured from the equatorial crossing, and a be the orbital arc 
measured also from the equatorial crossing. 



Figure ^ ^ 


Consider Figure 2.2 which illustrates the area element dA^ at the equator 



Figure 


2.2 
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It is obvious that dA^ is given by 

dA = d(|) de (2.1) 

o o o 

where r is the radius of the orbit. Let f denote the density of the satel- 

o 

lites at the equator. Then, the number dN of satellites contained in dA^ 
is given by 

dN = f dA (2.2) 

o o 

As these satellites move to latitude 0 and right ascension 4>> the corre- 
sponding area dA Is then given by 


dA = r COS0 d(f>d0 


and the density f is then obtained from 


(2.3) 


dN = f dA 


(2.4) 


Substitution of Equations (2.1) - (2.3). into (2.4) yields 


f d(j) d0 
o o o 

f = (2.5) 

COS0 d^d0 

However, from Figure 2.1, we obtain the following spherical trigonometric 
formula 

sin6 = sin i sin a (2.6) 


so that at latitude 6 we have 


cos0 d0 = sin i cosa do (2.7) 
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( 2 . 8 ) 


and at the equator we have 


d0 = sin i do 
o o 


Moreover, it is easily verified that we also have 


d<j) 

o 

-©- 

II 

(2.9) 

d a 

= d a 

(2.10) 


Substitution of Equations (2.7) - (2.10) into (2.5) yields 

f 


cos o 


( 2 . 11 ) 


which states that the density is inversely proportional to the cosine of 
the arc length. 

Next, we obtain the equatorial density f^ as follows: 


N = 


0 max C 2 ^ 

0 min o 
6max 


f r COS0 d^ de 


= 2Trr 


f cos6d0 


= 27rr 


= 27Tr 


o min 
. tt /2 

- 7T/2 
■ 7t/2 

- tt/ 2 


f sin i cos ado 


f sin i do 
o 


o 2 2 , . . 

2tt r f sin 1 

o 


( 2 . 12 ) 


in which Equations (2.7) and (2.11) have been used. 
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Substitution of Equation (2.12) into (2.11) yields 

N 


f = 

2 2 

2ir r sin 1 cos a 


(2.13) 


which expresses the density in terms of the total number, the radius, the 
inclination and the orbital arc. However, it is more convenient to obtain 
an expression in terms of latitude than orbital arc. This is accomplished 
as follows: Using the identity 


2 , . 2 
cos o - 1 - sin a 


(2.14) 


and also Equation (2.6), we obtain 


”2 2 

sin i cos a = /(sin i - sin 6) 


(2.15) 


so that Equation (2.13) becomes 


f = 


O 2 2 77 . 2. , 2„, 

/TT r /(sin 1 - sin 0) 


(2.16) 


2.2 Angular Separation 

The system of N satellites under discussion is considered to be uni- 
formly distributed in terms of equatorial crossing and also along the or- 
bital arc. It is easily verified that the angular separations between the 
satellites are given by 

(2.17) 

Ao = , (2.18) 
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SECTION 3 - REGION OF OBSERVATION 

3.1 Geocentric Conical Angle 

Consider a ground station G on the earth's surface. Let 3 denote 
the conical observation angle at the earth s surface, a the conical 
angle subtended at the earth's center, r^ the mean radius of the earth, 
h the satellite's altitude, and a the conical distance as illustrated 
in Figure 3.1. 



Figure 3.1 
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The geocentric conical angle a may be obtained as follows: For the triangle 

OGH, we have the sine formula 


a = 


(r + h) sin a 
e 

sin (ir - 6) 


(3.1) 


and the cosine formula 


a^ = r ^ + (r + h)^ - 2r (r + h) cos a 
e e e e 


(3.2) 


Substitution of Equation (3.1) into (3.2) and use of the identity 


.2 , 2 

sxn a = 1 - cos a 


(3.3) 


yield the following quadratic equation for cos a 


2 r 

cos a - 21 e 


sin 6 cosa + 


ir + h 
e 


^r + h 
e 


. 2 „ 

sin 6 - cos 6 


= 0 (3.4) 


whose solution is given by 


cosa 


2 T 

sin B + cosB 


+ h) 


1 -/^e 


r + hj 
e 


• 

sin 6 


(3.5) 


It may be verified that the physically acceptable solution is the one which 
yields the smaller angle a, i.e., the one with the positive sign in Equation 
(3.5). The other solution yields the larger angle a which results in a cone 

going into the earth, which is thus rejected. 

3.2 Boundaries of Intersection Region 

Consider Figure 3.2 which illustrates a ground station G and a ship 
S, and also the region of observation common to both of them. Let C be 
the central point of the great circular arc GS, and y the angle GOG. 
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Then, from spherical trignometry, the arc length 2y between G and S is 
given by 

cos 2 y - cos cos (3.6) 

and the inclination K of S with respect to G is given by 

0^ = sin K sin 2y (3.7) 

Next, consider Figure 3.4 which illustrates the boundaries R and L 
of the intersection region Q. It is to be noted that these boundaries 
are not arcs of great circles, but are arcs of small circles. 



Figure 3.4 

In order to obtain expressions for the boundaries R and L, it is conven- 
ient to consider the arc GS as the equator in an oblique coordinate 
system. First, consider the curve R. 
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Let 0 ' be the latitude and be the longitude of a point with respect 
to G. Then, from spherical trigonometry, the equation of the curve R is 
given by 


cosa = COS0' costjj" 


(3.8) 


However, if ii' denotes the longitude measured from C, then we have 


iP” = ij;' + Y 


(3.9) 


and Equation (3.8) becomes 


cosa = cos 0' cos (ij;' + y) 


(3.10) 


which is the equation for the boundary R in the oblique geographical system 
having C as the origin of latitude and longitude. Similarly, the equation 
for the boundary L is given by 


cosa = cos O' cos (tp' - y) 


(3.11) 


The 


points of intersection of the curves R and L are given by 


(ip> = 0, e' = e'max) and ((jj' =0, 0' = 0'min) where 


0 max 


-1 

cos / cos g 
cos Y 


( 3 . 12 ) 


min 


max 


(3.13) 
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3.3 Regular to Oblique Geographic Transformation 

Let ? = (x,y,z) denote the coordinates of a point in the regular 
geographic system, and r^= (x' , y' , z') denote the corresponding coordinates 
of the same point in the oblique geographic system. Figure 3.5 illustrates 
the angular rotations to accomplish the necessary coordinate transformation. 



Figure 3.5 


Let A denote the transformation matrix from r to r', i.e., 

r ' = Ar (3.14) 

Then, it is well-known that A is given by 



(3.15) 
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where the symbols s and c respectively denote the sine and cosine func 
tions of the argument which appears as the subscript. Next, it is also 
noted that r and r' may be respectively expressed in terms of their 
latitude and longitude as follows: 


X 

y 

z 





(3.16) 


(3.17) 


Thus, Equations (3.14) - (3.17) may now be used to express the oblique 
latitude and longitude in terms of the regular ones. The final results are 
given by 



s c.s , + c s. 

< 0 K 0 


(3.18) 


/ , \ p CiS , + s s. 

tan (ip + y) = k .9 k 9 


3.4 Oblique to Regular Geographic Transformation 

Next, to obtain the regular latitude and longitude in terms of the 
oblique ones, we proceed. as follows: We note that 

r=A^r' (3.20) 

T 

where A denotes the transpose matrix of A. 
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Then, proceeding as before but now using Equations. (3.20) and (3.15) - 
(3.17), we obtain 




'e' 


+ c 


' 0 ’ 


+ c 


^ 0 ‘ 


(3.21) 


tan tjj = 


+ Y) 


’e' 


S' + y) 


(3.22) 
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SECTION 4 


TIME AVERAGE POPULAllON 


4.1 Exact Folmulation 


Let N denote the number of satellites (time average population) 
within the common domain of observation It is obviously given by 


= 


4 


f dn 


(4.1) 


where the density f is given by Equation (2.16) and the element of area 
dn is given by 


dfi = r COS0 dtp d6 


(4.2) 


It appears that the above integral may be trivially expressed in terms of 
the regular latitude 0 and longitude tp as follows: 


0 max 



0 rain 


^( 6 ) 




N COS0 dij) d0 

2 2 2 
2Tf'^ y~ (sin i - sin 0) 


(4.3) 


where ip (0) denotes the expression obtained by solving for tp in terms of 0 
R 

using the equation for the R curve given by Equations (3.10), (3.18), and 
(3.19), and similarly for iri terms of the L curve. Not only is this 

process difficult, but it is noted that the integral on the RHS of Equa- 
tion (4.3) may not even be valid or, worse yet, amenable to numerical eval 
uation even in principle. This point becomes evident by combining Figures 
3.3 and 3.4. It is possible that the location of the ship S with respect 
to the ground station G can give rise to the case where, in performing the 
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integration with respect to ijj, the process does not take place from the L 
curve to the R curve and, furthermore, in performing the integration with 
respect to 6, the process also does not take place from 0 min to 0 max . 
This difficulty can be circumvented by writing the element of area df^ 
as follows 


dfi 


2 

r 


cos0 ' 


di|^' d0' 


( 4 . 4 ) 


so that the integral becomes 


0 ' max 



0 'min 





N COS0' d\|;' dO' 

2it^/~ (sin^ i - sin^0) 




r(6-) 

L(e').’ 



cosQ ' dif< ' de ' 


+ Vk^ 0' "ip- ^ Ve') 


■] 


(4.5) 


in which Equation (3.21) has been used. It is to be noted the 4^' integra- 
tion will always proceed from the L curve to the R curve, and the 6' 
integration will always proceed from 6' min to 0' max. 

4 . 2 Approximate Formulation 

An approximate formulation may be obtained by going back to the original 
Equation (4.1) which may be used to yield the following 



f n 

ave 


( 4 . 6 ) 
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where 


■# 


n = (ID r" COS0' dijj' d0' 

n 


= 4r 


rB max (q') 

Jo Jo 


COS0 ' dii^ ' d0 ' 


= 4r 


/ 


0 'max 


COS0 ' 


cos 


cosg 1 
COS0 ' i 


- Y 


d0' (4.7) 


in which Equations (3.10) and (3.12) have been used. This integral may be 
evaluated numerically once the relative position of the ship S is speci- 
fied. The average value of f to be used may be obtained by averaging the 
4 values at the mid-points on the axes of symmetry of 0. These, in turn, 
may be obtained by averaging the values at the center C and those at the 

extremities P illustrated in Figure 3.4. Thus, we may write 
1 


f 

ave 


1 

8 


f(P^) + f(P2) + f(P3) + f(P^) + 4f (G) 


(4.8) 
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SECTION 5 


AVERAGE POPULATION TIMES 


5.1 Exact Formulation 

First consider Figure 2.1, for which we may write the following spher 
ical trigonometric formulas 

sin 0 
cos a 

where the relevant quantities have already been previously defined in 
Section 2. Next, consider Figure 5.1 which illustrates the ground station 
at G, when the satellite crosses the equator at N^. Subsequently, when 
the satellite has moved to latitude 0, the rotation of the earth has taken 
the ground station to the point G. 



= sin 1 sin a 


= cos 0 cos <j> 


(5.1) 

(5.2) 


Figure 5.1 

Then, it is obvious that the following relation holds for both direct 
(i < tt/2) and retrogade (i > tt/2) orbits 

y + $ = V + ujo (5.3) 

where 

u) = — (5.4) 

P 

e 
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y “ longitude of satellite crossing measured from ground station 
<|) = right ascension of satellite measured from equatorial crossing 
ip = longitude of satellite measured from ground station 
o = orbital arc of satellite measured from equatorial crossing 
0 ) = ratio of satellite orbital period P to earth rotational period P 


Substitution of Equation (5.3) into (5.2) yields 

COSO = COS0 cos (l(j + 0)0 - p) 


(5.5) 


Equations (5.1) and (5.5) express the latitude and longitude in terms 
of the orbital arc. Symbolically, we may write 

6 = 0 (o; i) (5.6) 


ip = ip (a; i, p) 


(5.7) 


In turn, these equations may be substituted into Equations (3.18) and 
(3.19) to yield expressions for the oblique latitude 9' and longitude 
Ip' in terms of orbital arc o. Thus, we have 


0 ' 

<P' 


(0, 

i(); K, 

y) 

= 0' 

(o; i, y, 

K. y) 

(5.8) 

(6, 

\P; K, 

y) 

= >1'' 

(o; i, M, 

<, y) 

(5.9) 


which constitute 2 equations in the 3 unknowns 9', ip' and a. If we wish 
to determine the point of intersection with the R curve bounding one side 
of the common region of observation H, we also have Equation (3.10) which 
is . 

cos J. = cos 0' cos (i|)' + y) (5.10) 


Substitution of Equations (5.8) and (5.9) into (5.10) yields a com 
plicated equation for a which may then be solved numerically to obtain 
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the value corresponding to the intersection point. Next, to 

obtain the point of intersection with the L curve, we have Equation 
(3.11) which is really Equation (S.IO’) with Y replaced by -y . Thus, 
the same process may be used to obtain the value o=o corresponding 
to the intersection point. Thus, the population time T of the satellite 
within the region is exactly given by 


Let denote the value of y which corresponds to the orbit passing 
through the central point C. The above process is first performed with 
a value h " y(, + where is a random number in the range 0< A(j) < A<p 
where A(j) is given by Equation (2.17) which is ° 


(5.12) 


The process is then repeated with values (y + nAc))) where n = ±1, ±2,... 
until no more orbits intersect the region . After this, the entire 
above process is then repeated with other random values of Ac)) . The 

average population times are then obtained by averaging the re^sults of 
all these processes. 


5-2 Approximate Formulatio n 

Consider Figure 5.2 which illustrates the spherical triangle formed 
by the equator, the meridian and the arc length of the central point C 

measured from the ground station G. This spherical triangle is fixed on 
the rotating earth. 



q; 


Figure 5.2 
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Then, we have the following spherical trigonometric formulas 


sin0 

= sinK 

siny 

(5. 

13) 

cosy 

= COS0 

cost); 

(5. 

.14) 

siny 

sinX = 

sinij; 

(5^ 

.15) 


which may be used to compute the latitude and longitude of C and also 
the angle X the arc GC makes with the meridian through C. 

Next, consider Figure 5.3 which illustrates the spherical triangle 
formed by the equator, the meridian and the orbital arc of a satellite 
just passing through the point C. This spherical triangle is fixed on the 
celestial sphere, which is inertial. 


c 



(!) 


Figure 5.3 


Then, we have the following spherical trigonometric formulas 


sine 

= sin i 

sino 

Goscr 

= COS0 

cosc 

sine 

sine = 

s inc 


(5.16) 

(5.17) 

(5.18) 


which may be used to compute the orbital arc, the right ascension and 
also the angle ^ the arc NC makes with the meridian through C. 
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However, because of the earth's rotation, the satellite's ground 
track does not really make an angle ^ with the meridian through C. 
Rather, it is deflected through an angle 5 which is, in general, 
given by 

tan ? = (5.19) 

1 - o)cos i 


where co is defined by Equation (5. A). (It may be verified that this de- 
flection causes direct orbits to be more inclined and retrograde orbits 
to be less inclined as viewed by their ground tracks.) Thus, the angle 
between the satellite ground track and the meridian at point C is given 
by ( 5 -^), as shown in Figure 5.3. Next, consider Figure 5.4 which illu- 
strates the inclination r| of the orbital arc with the oblique equator 



Figure 5.4 

Hence, it is seen that we have 


n = A - C + ^ for ascending orbits 


n = A4-^-5-'rr for descending orbits 


(5.20) 


It may also be verified that these equations are algebraically valid for 
both direct and retrograde orbits. 
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Next, consider Figure 5.5 which illustrates the case of a satellite 
just passing through the point D whi^:h is displaced by Ay from the central 
point C. This corresponds to an orbit whose equator crossing is displaced 
by Acf) from the point N. 


D 



Figure 5.5 

Then, using spherical trigonometric formulas, it may be shown that Ay 
is related to A(() by the following equation 


tan i sin (t|j - (j) + A(j)) 

tan (y + Ay) = — (.5.^1; 

tan i cos k cos (\p - ip + Aip) - sin k 

Thus, by replacing y with (y+Ay ), Equations (5.13) - (5.20) may be used 

to compute the inclination n of the new orbital arc with the oblique equa- 
tor. It may be verified that Equation (5.21) is also algebraically valid 
for both the cases of i > K and i < K. Moreover , it is also valid for 
both direct and retrograde orbits. Furthermore, it is valid for arbitrary 
finite differences AcJ) and Ay , but considerable care must be exercised 
when taking the inverse tangent to obtain (y + Ay) in the correct quadrant 
corresponding to the increment Acf) . 
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Up to this point, no approximations have been made. It is now assumed 
that the satellite's ground track is an arc of a great circle lying within 
the region and making an angle r) with the oblique equator GS. Figure 
5.6 illustrates the cases of orbital arcs passing through the points C and D. 



Figure 5.6 


Now, it is possible to write the following two approximate relations for 
the orbit passing through the central point C 



where o’ is the arc length measured from the oblique equatorial crossing. 
These two equations are the crude analogs of Equations (5.8) and (5.9) of 
the exact case. If we wish to determine the point of intersection with 
the R curve, we also have Equation (3.10) which is 


cosa = cos 9' cos (ip' + y)- 


(5.24) 
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However, instead of substituting Equations (5.22) and (5.23) into (5.24) 
to yield a complicated equation for a', it turns out to be the case 
that an algebraic equation can be obtained involving sin 0'. This is 
accomplished as follows: From Equations (5.22) and (5.23), the following 

auxiliary equation is obtained 


sin ip' = 


tan 6 ' 
tan n 


(5.25) 


Equation (5.24) is then written as 


cosa = cos 0' 

(cos ]p' cosy 

sin ip' 

sin y) 

= cos o' 

cosy - cos 0 ' 

tan 0 ' 

siny 



tan n 


= cosy 

1l - sin^ 0’\ 

sin 0 ' 

siny 


\ sin^ n / 

tan n 



or equivalently 

2 2 

cos Y*^sin n - sin 6') = sinn cosa + cosn siny sin 0'. (5.26) 

By squaring both sides of this equation, it is obvious that a quadratic 
equation is obtained involving sin 6'. After much simplification, it may 
be shown that we have 


sin 0 ' 
sin n 

A little consideration will reveal that for the intersection point with 
the R curve, it is necessary to retain only the positive sign in the above 
equation. Thus, this expression corresponds to the value at 6 “ 9 r* 


/ — / .2 .2 . 2 - 

cosa siny cosn ± cosy / (sin a - sin y sin n ; (5.27) 

(1 - sin y sin n ) 
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However, from Equation (5.22), it is seen that the value o„ ' is 

K 

given by 


sin0 


R 


sinn sin o' 

R 


(5.28) 


Consequently, we have 


sin 


Next, to obtain the intersection point with the L curve, we have 
Equation (3;11) which is really Equation (5.24) with y replaced by -y . 
Thus, the same process may be used to obtain which can be shown to 
be given by retaining the negative sign in Equation (5.27). Consequently, 
we have the following result 

sin a ' 

Ju 


2 2 2 

= cosa siny cosn -cosyvTsin a - sin y sin'^n) (5.30) 


2 2 

(1 - sin Y sin n) 


, /7 . 2 .2 . 2 . 

cosa siny cosn + cosy /(sin a ~ sin y sin n) 


2 2 

(1 - sin y sin n) 


(5.29) 


which states that o' = - ' as expected (only for the case of the 

i-* K 

orbit passing through the central point C) . Thus, the population time 
T of the satellite within the region Q is approximately given by 

T = ^ (5.31) 

27T 

which is the crude analog of Equation (5.11). 

Next, to obtain the intersection point between the R curve and 
the orbit passing through the point D, a little consideration will re- 
veal that it suffices to replace y by (y + Ay) and also use the 
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corresponding value of ri and then repeat the process above for computing 

O ' as given by Equation (5.29). However, to obtain the intersection 

point between the L curve and the orbit passing through D, a little more 

caution is now necessary. It now suffices to replace Y by (-y + Ay) and 

also use the corresponding value of q and then repeat the process above 

for computing a ' but now retain the negative sign. This result yielding 
R 

the value of is no, longer trivially the negative of as for the 

special case of C. 

The above process is first performed with a random value in the 
range 0 < A<|> < A<j> where A(f) is given by Equation (2.17) which is 

A. 2’t (5.32) 

A(j) = — 

/n 

The process is then repeated with values (A(|)^ + nA(j)) where n = +1, +2,... 
until no more orbits intersect the region After this, the entire above 
process is then repeated with other random values of A(|)^. The average 
population times are then obtained by averaging the results of all these 
processes . 

Finally, it must be mentioned that in order to insure that the cor- 
rect segment of the R circle (see Figure 5.6) is identified to yield the 
desired intersection point as given by the general analog of Equation 

(5.29) , a little consideration will reveal that we must have n in the 

range -tt/ 2 < n 1 tt/ 2. Thus, if n is outside this range, we must accordingly 
add to or subtract it from n. Similarly, the same procedure applies to 
insure the identification of the correct segment of the L Circle to yield 
the desired intersection point as given by the general analog of Equation 

(5.30) . Furthermore, considerable thought will reveal that this assign- 
ment of the n range not only correctly gives the desired intersection 
points for orbits crossing the oblique equator within the observation re- 
gion n, but also for the case of equatorial crossings outside it for a 
range of A'v exceeding tt/2 measured from the central point C. The rea- 
sons for this are not apparent and, at first sight, it would seem that 
this assignment of n values outside the region Q leads to incorrect an- 
swers. But this is not so because of the manner in which the inverse 
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trigonometric functions are assigned their principal values. Thus, 
Equations (5.29) and (5.30) contain many subtle featpres in logic which 
automatically combine to yield, in mutual accord, the correct intersection 
points regardless of the equatorial crossing. In particular, additional 
consideration will reveal that it is only necessary to consider equatorial 
crossings such that the orbits intersect the oblique meridian through the 
central point C at an oblique latitude 9' not greater than 6* given by 



This corresponds to a range Ay* given by 


(5.33) 


so that 


Ay* 


. -1 
sin 


( tan 9* 

tan I n 



. -1 
sin 




(5.34) 


(5.35) 


or equivalently 


Ay* = min | tt/2, sin ^ /~/ 


2 2 
cos y - cos a 


(5.36) 


2 2 ^ 2 2 
cos y - cos a + cos a tan n 


It is not difficult to see that if an orbit intersects the oblique 
equator outside the range Ay* and also eventually intersects the observa- 
tion region Q, then this orbit would already have been counted as lying 
within the acceptable range on the other side of the central point. 
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SECTION 6 


RESULTS FOR ORBITING SATELLITES 


6-1 Average Population Time Computations 

Computations were performed, except for minor modifications, according 
to the method discussed in Section 5.2 to obtain the average population 
times for Class I and II satellites. The representative values of para- 


meters used are shown in Table 6.1. 

Table 6.1 

Quantity Class I Class II 

Period? (minutes) 100.9 717.9 

Inclination i (degrees) 74.0 63.9 

Altitude h (km) 800 - 20,178.5 

Number N 400 100 


The value of 3, the conical observation angle at the earth's surface, is 
taken to be 80° for both the ground station and the ship. The ground sta- 
tion is taken to be at the origin of longitude and latitude while the ship 
is taken to be at various values of longitude and lafitude 6^ only in the 
first quadrant. It may be verified that for locations of the ship in the 
other quadrants, the corresponding results may be obtained by taking mirror 
reflections about the primary axes. 

After the average population times t have been obtained, the results 
were divided by the characteristic time T defined by 

T = ^ (6.1) 

to obtain the number of satellites visible to both the ground station and 
the ship. (T is the time for a satellite to travel the intra-satellite 
distance Ao where Ac is given by Equation (2.18).,) The relevant results 
for Class I and II satellites are respectively summarized in Figures 6.1 
and 6.2, each of which was obtained by averaging the results using ^ given 
by Equation (5.i9) and those using ^ = 0. 
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Note: 


2 . 30% 

4 17% 


(0,30) 





(1) Numbers in the boxes denote 
the number of satellites vi- 
sible for the percent of time 
indicated. 

(2) Numbers below the boxes de- 
note the relative longitude 
and the absolute latitude of 
the ship. 
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97% 
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75% 

3 

65% 
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Figure -6.1 - Results for 100 Minute Orbiting Satellites 
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Note: (1) Numbers in the boxes denote 

the number of satellites vi- 
sible for the percent of time 
indicated. 

(2) Numbers below the boxes denote 
the relative longitude and the 
absolute latitude of the ship. 
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Figure 6.2 - Results for 12 Hour Orbiting Satellites 
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6 • 2 Time Average Population Computati ons 

For the special case of the ship at the origin of longitude and lati- 
tude, the time average population may be computed by Equation (4.3). 
Numerical integration yields a value of about 28.48% for N /n for Class II 
satellites. That is, on the average, 28.48 satellites are mutually visible 
to the ground station and the ship when they are together. 

As a comparison, it may be shown that the ratio of the area of common 
visibility SI to the area of the zonal belt A covered by the satellite 
orbits is given by 

^ [l 
A 

when the ground station and ship are together. Hence, for Class II satel- 
lites, we obtain a value of about 22.4% for fi/A. As expected, this value 
IS smaller than that for N^/N because the density f increases with latitude 
and hence contributes toward giving a higher value of in the numerical 
integration. 

The other comparison is made with the results displayed in the (0,0) 
box of Figure 6.2 which are seen to yield a smaller value than that for 

is also to be expected because the approximation made in Sec- 
tion 5.2 assumed that the satellite orbits are arcs of great circles within 
the region and hence yields a smaller value of the average population 
time T than that obtained by considering the actual satellite ground track. 
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DISTORTION-FREE MAPPING OF VISSR IMAGERY 
DATA FROM GEOSYNCHRONOUS SATELLITES 

F. K. Chan 

Scientific Analysts and Consultants, Inc, 
lillU Heathfield Road, Rockville, Md, 20853 


ABSTRACT 


Analysis has been performed for mapping VISSR imagery 
data so as to eliminate all geometrical distortions. The formula- 
tion is rigorous and includes all misalignment angles of the VISSR, 
the sun sensor and the instantaneous spin axis with the satellite's 
body axis. It also includes the effects due to the motion of the 
satellite's suborbital point. All the mapping equations for dis- 
tortion removal are reduced 'to simplest forms, and all the algorithms 
are optimized as much as possible. 

An approach is then formulated for implementing these 
algorithms for in-line operational use . It covers the computations 
involved in determining benchmarks, the interpolation methodology 
for filling in the points interspaced between benchmarks, and the 
correction procedure for computing the radiometric values at the 
center of the pixel in the distortion-free image. It is also concerned 
with the time requirements, data storage, and output data accuracy. 
With the present microprocessor technology, it is concluded that this 
in-line distortion removal is possible in real-time processing of 
infra-red but not visible VISSR imagery data* 


This work was supported by NOAA. Contract Nos; 01 -8-M01 -1 86U 
and NA-79KAC-00026 
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SECTION 1 - INTRODUCTION 

In the Visible and Infra-red Spin Scan Radiometer (VISSR) 
data obtained from the present geosynchronous satellites, distortions 
are observed in the images of the earth. As illustrated in Figure 1 .1 
which is exaggerated for clarity, these image defonnations appear 
as vertical compression and expansion of the image, non-vertical 
alignment of the North and South Poles, and multi -representation of 
some points or omission of other points. 



Figure 1.1 Exaggerated VISSR Image of Earth 
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These distortions may well be explained by considering 
Figure 1.2 which illustrates the same scan-lines; on the projection 
plane of the earth. Again, for clarity, these scan-lines are depicted 
to be non-parallel and unevenly spaced to a degree more so than the 
realistic cases. 



s 


Figure 1.2 Projection Plane Image of Earth 
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If one were to relate these two images, one would find, 
for example, that a triangijlar figure in the projection plane 
image becomes distorted into a curved figure in the VISSR image. 

This is illustrated in Figure 1.3 which is obtained by super- 
imposing a triangle on Figure 1.2 and then mapping it onto Figure 1.1, 

z 

5 

4 

S 

Projection Plane Image 





The general causes for these image deformations may be 
broadly classified as follows: 

1 . Orbit not circular and equatorial 

2. Spin axis not perpendicular to orbital plane 

3. Misalignment of the VISSR, the sun sensor, and the 
instantaneous spin axis with the satellite's body axis 

U. Biases due to varying sun size and varying sun elevation 
effect on threshold of the sun sensor triggering. 

To remove these distortions, it is necessary to include all the 
above factors in the formulation of the mapping equations. However, 
it is feasible to consider only the first three. The corresponding 
equations have been derived in Reference 1 in which it was convenient 
to introduce the following coordinate systems : 

The Inertial System : This is well-known and is defined such that 

the Xj-axis is in the direction of the vernal equinox, the Zj-axis 
is perpendicular to the equatorial plane (in the direction of the 
North Pole), and the yj-axis is given by y^. ■ z‘j x x^. 

The Body System : This system is defined such that the Zg-axis is 

along the longitudinal axis of the satellite, the Xg-axis is the 
intersection line between the VISSR stepping plane and the plane 
perpendicular to the Zg-axis, and the yg-axis is given by yg « Zg x ^ 
The VISSR System : This system is defined such that the x^-eucis is 

in the direction of the mid-scan, the z^-axis is perpendicular to the 
x^-axis and lies in the VISSR stepping plane, and the y^-axis is 

, , A /% 

given by y^ « Zy X x^. 

The Sun Sensor System : This system is defined such that the Xy-axis 
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The Auxiliary System ; In this system, illustrated in Figure 1«li, 
the unit base vectors are defined by the following equations: 



Figure 1 ,U - The Auxiliary System 
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The Normalized System ; Consider a system, referred to as the normalized 
system, as illustrated in Figure f.$. The origin R of this system is 
defined at a point"? on the earth's equatorial plane and fixed in the 
earth's rotating system. The satellite P at point r', however, is not 
necessarily on the earth's equatorial plane or fixed in the earth's 
rotating system. In this normalized system, the unit base vectors 


are defined by the following equations: 

A A 

* - A 

A A 

A A A 


ChO 

C('i) 

0 « 8 ) 
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In the present study, the results obtained in Reference 1 
are used to formulate algorithms for mapping the data coming out of 
the Synchronous Data Buffer (SDB) so as to obtain a distortion-free 
imagery. Moreover, this rectified imagery also has the desirable 
feature that it is referred to a normalized satellite position which 
is therefore the same for all imageries. Thus, if the distortion-free 
mapping is performed in-line during data processing, the transmitted 
VISSR data will provide a uniformly compatible data base for all 
users in their scientific work. Furthermore, it will also facilitate 
in the future development of a composite data base for different 
kinds of data obtained frcm various satellites. 

Mapping of the data may be further optimized the use of 
interpolation with the aid of appropriately chosen benchmarks. 

Section 2 deals with the computation of these benchmarks, while 
Section 3 covers the interpolation methodology for filling in the 
points interspaced between benchmarks. Section U is concerned with 
the correction procedure for ccmiputing the radiometric values at the 
center of the pixel in the distortion- free image. Section 3 discusses 
the time req\xirements , data storage and output data accuracy. 

Section 6 summarizes the results of this study. 
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SECTION 2 - BENCHMARK COMPUTATIONS 


This section deals with the computation of benchmark 
locations in the normalized distortion-free coordinate system. It 
discusses the relevant input parameters, computational equations, 
number of compuational operations, and the requisite partial deriva- 
tives. 

2.1 Input Parameters 

The relevant input parameters are listed below; 

N,j^ = scan number 

N, = sample number 

= mid-scan, number 

Mp = mid-sample number 
"is 

A, = scan angular width 
= sample angular width 
“ scan angular bias (line bias) 
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"r' = orbital position of Satellite 


Oi = right ascension of spin axis 


5 = declination of spin axis 


"t - normalized position of satellite 


a = earth's semi -major axis 
c = earth's semi -minor axis 

For convenience, two parameters dependent on tne above are defined 
as follows: 
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2.2 Computational Equations 

The following equations for computing benchmark coordinates 
have been extracted from Reference 1 . They have been simplified and 
are listed below in 'the proper sequence for usage. The exact definitions 
of cursory intermediate variables may be obtained from the original 
report. 
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I" J’"' are the coordinates in the normalized distortion-free 
system. At this stage, for the sake of greater accuracy in sub- 
sequent computations, it is preferable not to digitize them. 
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2.3 Computational Operations 

For each benchmark, it may be verified that the comput- 
ations in equations (2.1) - (2.28) may be achieved by performing: 

additions 
22 divisions 
11 ii multiplications 
ii divisions 

1 1 trigonometric function evaluations 
2 square root evaluations 

Assuming that the following times are required; 


Operation 

Time (micro-seconds)' 

Addition 


Subtraction 

1.5 

Multiplication 

6 

Division 

11 

Trigonometric function evaluation 

50 

Square root evaluation 

50 


it is seen that about 1U87 microseconds are required for each bench- 
mark computation. 
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2.ii Partial Derivatives 


For convenience, it is desirable to choose the set of 
benchmarks so that they form a rectangular grid in the (n^ > N ) 

U 'I 

-space as illustrated in Figure 2.1. 

1 

Mr* • s 


p • 


Figure 2.1 - Benchmarks in (N_ > N ) - Space 

X 


Then, it is obvious that the partial derivatives of 1“ and 
with respect to N> and N may be approximated by 

Ss n 


'hN. 
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91 


9/V, 





(i% - (r)p 
(r)i, - (r)p 


( 3 ^ 31 ) 


C3^B2) 


For each benchjnark, the four associated partial derivatives 
require 6 subtractions and h divisions. These operations consume 
about 53 microseconds. 


These partial derivatives are used later in the method of 
interpolation for mapping points interspaced betweeen the benchmarks. 
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SECTION 3 - INTERPOLATION COMPUTATIONS 

This section is concerned v/ith the mapping of points 
which do not coincide with the chosen set of benchmarks. It 
discusses the input data, interpolation methodology, and number of 
computational operations. 

3.1 Input Data 

The input data consists of the coordinates (T^, and 
their four associated partial derivatives for each benchmark. This 
information has already been obtained in Section 2. 


3,2 Interpolation Methodology 

Suppose there are m interspaced points between the hori- 
zontal benchmarks, and n interspaced points between the vertical 
benchmarks. Figure ).1 illustrates a basic unit comprising bench 
marks (denoted by solid circles) and interspaced points (denoted by 
open circles). Thus, there are (m + 1)(n + 1) points altogether 
in a basic unit. 
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Figure 3*^ " Benchmarks for Interpolation 


It is noted that the I" coordinate of the point T directly 
below P is given by 

- -I- ^ i(^\X 

and AH^ ° , 


(Ssi) 

^ 3 . 2 ) 
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Therefore, it follows that 


{r\ = (i% ] 


^3*4) 


, P 


A similar equation holds for any point and the point directly 
below it. Hence, the n interspaced points between P tod R may be 
mapped by the following iterative algorithm: 


Let 


if = 


Cb>s-} 


If - (A 




Then, perform the following computations 


I 




/ 


J 




I 4 
I 


T* 


mu 


1 

* 


[fUl) 

LLdfl^ 




0 A A, 




Similarly, it is noted that the I" coordinate of the 
point U directly to the right of P is given by 




(2^9) 
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But 


and 




S5 


A Ny = •/• / 


(i*lo) 




Therefore, it follows that 


^ mxx 


( 3 -/ 2 ) 


Again, a similar equation holds for any point and the point directly 
to the right of it. Hence, the m(n + 1) interspaced points between 
the columns PR and QS may be mappea by the following iterative 
algorithm: 


Let 


0,1 


r 

'^0,1 


1 -= 0 , 1 , 1 , . . . n 




(3. /a") 




Then, perform the following computations 



1 

ff¥) 1 ] 

1 

(■k+i ), t 



1 

L J ^ ^ 



t n 


J* V- 

j 

[fill] 7 J 

\ 

( 3 ^ 0 ,') 


After these computations, it is now permissible to digitize I” and J 
to obtain the rounded integer values I and J. 
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3.3 Computation Operations 


For each interspaced point, it is seen that the comput- 
ations in equations (3.7) and (3.8) or those in equations (3.15) 
and (3.16) require 2 "additions. Assuming a time of 1.3 microse- 
conds for each operation^ "therefore about 3 microseconds are required 
to map each point by interpolation. Allowance is also to be make for 
converting tv/o real numbers to integer values for each point. This 
will probably increase the time requirement by a factor of 2 so 
that about 6 microseconds are required for each point. 
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SECTION U - RADIOMETRIC COMPUTATIONS 


This sections considers the methodology of correcting the 
radiometric values so as to reflect a more realistic value at the 
center of the pixel in the distortion-free image. It discusses the 
input data, correctidn methodology, and number of computational 
operations. 

ii.1 Input Data 

For each point, the input data consists of the follov7ing 
and may be obtained either from the SDB output data stream or has 
already been obtained in Section 3: 

N^ = scan number 

a 


N- = sample number 


R(lt , N ) = radiometric value of N^ th scan and N_ th 

’I ’I 

sample 



) = radiometric value of th scan and 
(N- - 1 )th sample 


R(N _ , N -1) = raidometric value of (N - 1 )th scan 

>$ T 1 

and Nj, th sample 
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I = interpolated horizontal coordinate, -of 
in distortion-free system 




interpolated vertical coordinate of 
in dis'tortion-free system 



) 


-X- 

I = rounded integer value of I 


J = rounded integer value of 


I 4 . 2 Correction Methodology 


value R(I 


The radiometric value R(I, J) may be obtained from the 
, J") by using the Taylor's series expansion 


(I, 
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/dR N 

The oartial derivatives ( ^TF / *and 
mated by- 



may be approxi- 






- [ThX 




The four remaining partial derivatives on the RH3 of equations [h»2) 
and (U.3) may be obtained as follows: Let A ggid B be matrices defined 

by 



7-25 



Then, from the theoi*y of matheraaticaX transformations, we have 

A = E>~' 


which explicitly yields 



The partial derivatives appearing on the RHS of equations (U.9) - 
(U.13; rnay be obtained from equations (2.29) - (.2.32), valid for 
basic unit defined by benchmarks. 


(4-i) 

(’f'/o) 

(i'll) 

(^• 13 ) 

( 4 -/ 3 ) 
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U. 3 Computational Operations 


if’or each point, it is seen that the computations in 
equations (U. i / - (U.5) require h additions, 6 subtractions and 
6 multiplications. Assuming a time of 1.3 microseconds for each 
addition or subtraction , and 6 microseconds for each multiplication, 
therefore about 31 microseconds are required to correct the radio- 
metric value for each point. 
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SECTION $ - DISCUSSION 

This section discusses relevant topics such as benchmark 
spacing, time requirements for mapping IR data, real and non-real time 
computations, data storage and buffering, input and output data ac- 
curacy, and computational accuracy requirements. 

5. 1 Benchmark Spacing 

The IR samples have angular widths of about 0.01° x 0.005° 
at the satellite position. This corresponds to a resolution of about 
Ii X 2 miles at the sub-satellite position on the earth's surface. In 
general, this resolution and the non-linearity of the mapping equations 
determine the requisite spacing of the benchmarks to be used for inter- 
polation. The analytical approach to obtain this spacing involves 
comparatively complex mathematical aneilysis. Alternatively, it is 
also possible to obtain this value by actually performing the mapping 
munerically. At this stage, it is felt that the interpolation require- 
ments can be met by choosing the IR benchmarks to be spaced ^0 samples 
horizontally and 2$ samples vertically. That is, it probably suffices 
to choose m = 50 and n = 25 in Section 3. In the full IR imagery, 
there are 1822 scans each containing 3822 samples. Consequently, about 
5,600 benchmarks will be required. 
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5.2 Time Requirements 

In Section 2, the computation of each benchmark and its 
associated partial derivatives requires about ^$h0 microseconds. 

Hence, a set of ^600 benchmarks requires about 8,62h,000 microseconds 
— 8.6 seconds. 

In Section 3> the mapping of each sample by interpolation 
requires about 6 microseconds. Therefore, an IR imagery of about 
7 X 10^ samples requires about h2 seconds. However, if the entire 
IR imagery is not to be mapped, then cropping out the edges will 
probably reduce time by a factor of 2/3 to yield a requirement of about 
28 seconds. 


In Section U, the correction of radiometric value at 

the center of the pixel in the distortion-free image requires about 

5l microseconds for each sample. Therefore, an IR imagery of about 

6 

7x10 samples requires about 357 seconds. Cropping will probably 
reduce this to about 238 seconds. 

Consequently, about 399 seconds or 6.7 minutes will be 
needed to map the entire IR imagery cmprising of coordinates and 
radiometric values of the samples. This time requirement drops to 
about U.U minutes if cropping is introduced. 


7-29 



If it is desired to map* the visible imagery containing 
1/2 X 1/2 mile samples, then the above times are increased by a 
factor of 32. Therefore, about 2lU minutes will be required to map 
the entire imagery comprising of coordinates and radiometric values. 
If the edges are cropped out, then about II4.3 minutes will be needed. 


5.3 Real and Non-Real Time Computations 

From the discussion above, it is seen that it is possible 
to perform all the mapping computations in real-time in the case of 
IR imagery, and not possible in the case of visible imagery. However, 
in the latter case, the crucial point is whether the radiometric 
corrections are really necessary. If not, then the time requirements 
drops to 22. U minutes for the entire imagery, and lU.9 minutes for 
the cropped imagery. Consequently, visible data-mapping becomes feas- 
ible in real-time. 

Because the benchmark computation time is so small, it 
is desirable to perform the benchmark computations in real-time so 
that the relevant parameters may be easily extracted in-line from 
the data-stream coming out of the SDB. 
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5.U Bata Storage and Buffering 

Because the imagery obtained from the SDB output data is 
distorted, it is necessary to store this data in a buffer before the 
distortion-free mapping can be performed. The buffer size may be 
estimated by allovjing for a maximum 3° offset in the spin axis. Since 
the satellite is about 6.6 earth radii away, it may be verified that 
about 100 IR scan-lines (382,200 samples) to be buffered at a time. 
This will be sufficient to output a horizontal distortion-free line 
from end to end. In the case of visible data, the corresponding 
buffer will contain about 800 visible scan-lines (12,230,itOO samples). 
If a realistic situation, the above numbers will probably be reduced 
by a factor of 3* 

5.5 Input and Output Data Accuracy 

The data coming out of the SDB will be used as input 
into the distortion-free mapping software system- The acciiracy of 
this data may be roughly classified as perfect, normal or bad. 

Perfect data corresponds to data having errors of less 
than one pixel (i. e. ♦ 2 km at the subsatellite point for IR data). 

The error in the output data from the distortion-free mapping is 
therefore determined by the pixel resolution of the benchmarks, the 
interpolation accuracy/ of interspaced points, and the correction accur- 
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acy of the radiometric values. Because the second and third factors 
depend on the benchmark-spacing, which in turn de^iends on the pixel 
resolution, therefore it is estimated that the error bound of the 
output data is about one pixel (i.e., + U km for IR data). 

Normal data corresponds to data having errors of about 
one or two pixels. The mapping error is determined by the benchmark 
accuracy corresponding to normal input error, the interpolation 
accuracy of interspaced points, and the correction accuracy of the 
radiometric values. In this case, the error bound of the output 
data is about two pixels. 

Bad data corresponds to data having errors of about U or 
more pixels. The mapping error is determined mainly by the benchmark 
accuracy corresponding to these bad input errors. In this case, the 
error of the output data is probably about ^ or more pixels. 

5.6 Computational Accuracy Requirements 

It is desirable to investigate into the use of l6-bit 
words in the distortion-free computations. 

because of the complexity of the benchmark computations 
in equations (2.1) - (2.28), it is quite evident that sufficient 
accuracy will not be obtained by performing single-precision 
computations using 1 6-bit words. 
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However, for the partial deri”atives comnutations in 
equations ( 2 . 29 ^ - ( 2 . 32 ;, the interpolation of interspaced points 
computaions in equations (3.7; - (3.h) and (3.1?) - (3.16), and the 
radiometric correction computations in equations (U.1) - kU-?), 
it is possible oo achieve the desired accuracy using single-precision 
computations involving 1 6-bit words. In this case, perhaps the best 
way to represent real numbers is as follows: 

I bit for sign of nimiber 

1 1 

I I bits for range of number (2 - 1 = 20h7 ) 

1 bit for sign of exponent 

3 

3 bits for range of exponent (2-1=7) 

An alternative choice is as follows: 

1 bit for sign of number 

12 

12 bits for range of number (2 - 1 = U09? ) 

1 bit for sign of exponent 

2 

2 bits for range of exponent (2 -1=3) 

This second choice may not be as desirable because of the small 
range of the exponent. 
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SECTION 6 - CONCLUSION 


From the preceding discussion, it is seen that it is 
possible to map in real-time the entire IR imagery comprising of 
coordinates and radiometric values of the samples. However, it is 
possible to map in real-time the entire visible imagery comprising of 
only the coordinates of the samples. This conclusion is based heavily 
on the assumption that it takes 1.5 microseconds for each addition or 
subtraction, and 6 microseconds for each multiplication. 

The follox«.ng table summarizes the time requirements for 
IR and visible imagery mapping: 



Entire IR 

Cropped IR 

Entire VIS 

Cropped VIS 

Benchmarks 

8.6 sec. 

5.7 sec. 

U.6 min. 

3.1 min. 

Samole 

Coordinates 

U2 sec. 

28 sec. 

22. U min. 

1h.9 min. 

Radiometric 

Values 

357 sec. 

238 sec. 

I 90 .U inin. 

126.9 min. 

Total 

6.8 min. 

U.5 min. 

217 . h min. 

IUR .9 min. 
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The following table summarizes the expected accuracy of 


distortion-free mapping 

(DFM) algorithms: 

Input Data 

Output Data 

from SDB 

from DFi'^ 

Perfect 

1 pixel 

Normal 

2 pixels 

Bad 

5 or more pixels 


REFERENCES 

(1) Chan, F. K., "Distortion-Free Mapping of VISSR Imagery Data 
from Geosynchronous Satellites", Scientific Analysts and 
Consultants, Inc. Report. (1978). 


7-35 




Computational Aspects of Geometric 
Correction Data Generation in the 
Landsat-D Imagery Processing.* 

I . Levine 

General Electric, Space Division 
4701 Forbes Blvd. Lanham, MD 20706 

ABSTRACT 

A method is presented for systematic and geodetic correction data calculation. 

It is based on presentation of image distortions as a sum of nominal distortions 
and linear effects caused by variation of the spacecraft position and attitude 
variables from their nominals. The method may be used for both MSS and TM image 
data and it is incorporated into the processing by means of mostly offline 
calculations. Modeling shows that the maximal errors of the method are of the 
order of 5m at the worst point in a frame; the standard deviations of the average 
errors less than .8m. 

INTRODUCTION 

The geometric correction of the Landsat-type imagery typically proceieds in 
two steps. The first, called the Systematic Correction, removes internal distortions 
Imported in the raw image data by the sensor mechanism, spacecraft motion, inaccurate 
sensor pointing, earth's rotation, etc. These partly corrected Images still contain 
distortions due to uncertainties in spacecraft position and orientation. The second 
step. Geodetic Correction, removes these residual distortions using refined values 
of the attitude and ephemerls estimates. The refined attitude and ephemeris are 
obtained by filtering of image dislocations at Control Points. 

Application of the geometric correction requires the generation of the Correction 
Data - Systematic (SCD) or Geodetic (GCD) , depending upon the processing step. 

This data is developed on a rectangular grid in input (pixel, scan line) coordinates 
and express the relationship between the input and output map coordinates, within 
a standard World Reference System (WRS) frame. 

The central part of the SCD/GCD generation Is the computation of the coordinates of 
the intersection of the sensor's llne-of-slght vector, with the Earth's surface 
(lookpoint coordinates). The lookpolnt coordinates must then be converted to 
geodetic coordinates followed by mapping into user's map projection. There are two 
user's map projections: Space Oblique Mercator (SOM) and either Universe Transverse 

Mercator (UTM) or Polar Stereographic (PS) . 


* Work performed under National Aeronautics and Space Administration 
Contract No. NAS 5-25300. 


8-1 



Finally, the data, computed for integer values of pixels and lines, is inter- 
polated to integer values of output map coordinates. The grid spacing is chosen 
so that the data, together with properly defined interpolation techniques, represent 
the output coordinates to the desired precision everywhere in the frame. 

It should be noted that all the calculations are performed twice at each grid 
point, once for each SCD and GCD. They consume a significant amount of the processing 
time, which needs to be minimized. At the same time, there are no essential 
differences between SCD and GCD. Both establish a pointwise transformation, which 
may be written generically as 


^m ~ F(pixel,line,p) , 


where ^m2^ coordinates of a grid point and p is a vector of 

variables characterizing the spacecraft motion, attitude pointing, sensor's 
mechanism, etc. 

Letting p = p’^ + 6, the sum of nominal values of the variables and the deviation 
from the nominals, in the first approximation 


X = X 
m m 


+ y 6, 


( 1 ) 


where X^ are the nominal map coordinates and y is the matrix of the partial 


derivatives (PD) 



Thus, SCD and GCD may be represented as a sum of the nominal correction data and 
pointwise adjustments. 


This has significant advantages: 


1) It provides -a uniform approach to the SCD and GCD computations, considering 
each as one transformation, and 

2) Because the nominal spacecraft motion is known for every WRS frame, the nominal 
coordinates and the partial derivatives may be computed and stored in a Data Base. 


The Implementation of such an approach depends a great deal on both the choice 
of an output map projection and 6. An analytic form of mapping not only has to 
allow derivation of the coefficients y, but it should also afford rapid and precise 
online inversion to geodetic coordinates, from which the final map projection can 
be generated. In addition, it is desirable to have the nominal coordinates and the 
partial derivatives, as far as possible, insensitive to global position of the 
frame. Thus, although out technique may be applied to most standard map projections 
(such as UTM or PS), a special intermediate projection. Local Space Oblique Mercator 
(LSOM) , has been employed. The LSOM is the Mercator projection for the sphere, with, 
local 'equator' along the nominal spacecraft inertial velocity at the frame center. 
In that projection X and y are longltude-undependent and thus, can be stored only 
for one path. A natural choice of variables 6 is the along-track, cross-track and 
radial deviations in spacecraft position, together with deviations in the attitude 
angles. The nominal spacecraft motion within a frame is assumed to be in a perfect 
circular orbit passing through the frame center. 
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NOTATIONS 


ap (®pl 5 ^p2 » • ' • ’ 

II^pII 


— vector (Kxl matrix) 


the Euclidean norm of a 


A 


transposed matrlcs A 


X 


spacecraft position vector in earth- 
centered earth-fixed coordinates 


X 


spacecraft position vector in nominal 
spacecraft coordinates 


X = 
m 




- output map coordinates 


pointing vector in body coordinates 


pointing vector in local vertical 
spacecraft coordinates 


pointing vector in earth-centered fixed 
coordinates 


u 


coordinates of a lookpoint on the earth 
surface 


distance from spacecraft to earth 
lookpoint 


V 


u • |u|| ^ - normalized lookpoint vector 


R 


- local earth radius at WRS frame center 


earth rotation rate 


^e’ ^e 


ej^ ®2 ^e ^e^e 


earth equatorial and polar radii 

-1 
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deviation in the pitch 



deviation in the roll 

deviation in the yaw 

along track angular deviation 

cross track angular deviation 

relative departure in the radial direction 

time deviation 

(2x1) matrix of the partial derivatives 

of X with respect to 6, . 
m k 


P = (PpP 



a 

s 


X 


2 > • ■ • »P y) 

the 'equivalent' pitch and roll 
spacecraft orbital rate 
geodetic longitude 
geocentric latitude 


f - geodetic latitude 


ROTj^('F) = 


1 

0 

0 


0 

cost|) 

sinij; 


0 

-sini|) 

cosi|; 


ROT^Cijj) = 


ROT^Cij;) = 


cosif/ 

0 

sini|) 

0 

1 

0 

-slni|) 

0 

cosif) 

cosi(; 

_ 

-sinif) 0 

sinif) 


cosif) 0 

0 


0 1 
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J- - 0 ■ 

0 0 0 
^10 0 -1 
0 10 

0 0 1 
0 0 0 
-1 0 0 

0 -1 0 
T = 1 0 0 

■^ 0 0 0 


I - the three dimensional identity matrix 

The upper index n indicates the nominal value of a vector. 

T ^ - active scan time 
act 

T - mirror turnaround time 

round 


THE NOMINAL SCO 
Coordinate Transformations 

The local (instantaneous) spacecraft coordinates are described in terms of 
the unit vectors where 5^ points towards the Earth center, the 

vector is along the orbital angular momentum, and '$ 2 - lo ^ "Ci is roughly along 
the velocity direction. The local spacecraft coordinates at the WRS center is 
called the nominal anacecraft coordinates. The ma_trix A transforms a vector X in 
earth-centered inertial coordinates to the vector X in nominal spacecraft ° 
coordinates: ® 

X = AX (2) 

so ' 

The inertial to earth-centered earth-fixed coordinate transformation is 
defined as 

X = e \, (3) 

where E = E • R0T_(f2t). 

o 3 

The matrix gives the time-independent component of the transformation, 

ROT^(J^t) describes the rotation about the earth axis at the rate fii. We assume that 
t = 0 at the frame center. The corresponding nominal spacecraft to earth fixed 
coordinate transformation may be written as 

— T T— T — 

X = E A X = P X (4) 

s s 


where P = AE. 
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The unit vector g’^, given in body (sensor) coordinates, is transformed to local 
spacecraft coordinates as 

g = ROT^C-Oy) • ROT^C-e^) • ROT(-Qp) (5) 

where 0,0,0 are the yaw, roll, and pitch. 

y> r’ p > 5 r 

In the nominal spacecraft coordinates, g^may be expressed as 
is = ROTj^(Y)g 

where Y is the angle in the orbit plane between the spacecraft and the f^ame center. 
In the nominal spacecraft motion cos Y = - XsS/ l^sll » sinY = ^3 2^1 ’ 

and thus, the matrix ROT^^C Y) is known completely. 


A vector X = (X , , X „) in LSOM coordinates is defined as 
m ml’ m2 


X 


ml 


R , 1 + sin 

-7T In -Y : — 

2 1 - sin 


(7) 


X „ = Ra 

m2 


where R is the earth local radius at the frame center. The local polar angles, 
a and g , are given by 


= sing 

W2 = cosg • sina (8) 

= -cosg • cosa 

where _ _ 

W = II u II Au (9) 

and u = (u^ ,u ,u ) are earth fixed coordinates of the corresponding point on the 
ground. 
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Generation of the nominal SCD 

The nominal coordinates, Xm, are computed on a grid, consisting of 2xx\+l 
fictitious scan lines, each line containing 2n2+l points. Because the TM scans in 
two directions, it requires two sets of the nominal coordinates, for forward and 
backward scans. The computations may be fulfilled in the following order. 


1. Generate the time of (l,j) point 


t , 


li 2 n 


. q . C . e . ne . (i-nj^-1) + ■ (j-n^-1) + AT 


^n2 


Here 


T = (T ^ + T .) • (N - K.) 

scene act round scan 1 


where T is the active scan time, T a turnaround time, N is the 

actual number of scans, and = 1 for°MSS and 2 for TM. scan 

The parameter T = T for backward scans of the TM and zero otherwise. 

scene 

2. Generate 2n2+l unit line-of-sight vectors g’^ in body coordinates. An actual 
mirror velocity profile, together with constant sensor's misalignments may be employed 


3. Compute the spacecraft position vector X and the matrix P at t . 

s ij 

4. Compute g^ according to (6) . 

5. Transform X^ and g^ in earth fixed coordinates, obtaining the vectors X and f, 
respectively. 

6. Determine the lookpolnt coordinates, u = (u^,U 2 ,U 2 ), and h from the equations 

u = X + hf (10) 

(uj + u^) a-/ + . 1 (11) 

7. Transform u into LSOM coordinates using (9), (8), and (7). 

It is convenient to have all distances in units of the nominal orbit radius. 


THE PARTIAL DERIVATIVES 

Position and Pointing Vectors 
~*^ri 

Let X^ be a nominal spacecraft position vector at time t , 5 and 6 be the 
angular aiong-track and cross-track deviations in spacecraft position, and 5 be 
a relative deviation in the radial direc_tlon. Then the actual spacecraft position 
vector, X^, may be obtained by rotating X'^ through 6, and 6 . This is followed by 

stretching according to the ratio 1 - 5^:® ^ ^ 

6 

Xg = R0T2(<S3) • ROT^(d^) • X^(l-5g) 
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Similarly, if ^3 deviations in the pitch, roll, and yaw, the actual 

(unit) pointing vector should be written as 

is= ROT 2 ( 65 ) • R 0 T^( 6 ^)i" , 
where i" is defined by (5) and ( 6 ) as 

= ROT^(y) R0T2(-52)R0T2(-<52) • R0T^(-5p i’^ 

Let t = t + 6 -, . Remembering that P = P(t ) = AE • ROT„(J^t ) , we can write 

Q / O O J AJ 

T . 

P at txme t as 

p'^(t) = ROT^ (-fi(t^+ 6 ^))E^'^A'^ = ROT^C-S^d^)’ 

•ROT„(-S^t )e'^a'^ = R0T„(-fi6 )P^ 

3 o o 3 / 


and the actual position and pointing vectors in earth fixed coordinates as 
X = R0T2(-Q67)P^R0T2(§5)R0T^( d^)^(l-<Sg) 
f = R0T2(-fi6^)p\0T2(d3)R0T^(d^)R0T^(Y) , 

ROT2(-d3) • R 0 T 2 (- 62 ) ‘ ROT^C-e^)!’^ 

The linear terms of the Taylor series expansions of X and f in the vicinity of 
6^ = 0 (i = 1,2,. . . ,7) give 

X = pT (I + 1 ^ 6 ^+ 1265-6 

f = p'^ j^I+Tj^ROTj^(Y)d^+T2ROT^(Y)d3-ROT^(Y) • T^6j,- 

- 5^ • R0T^(y)6 J 


Here we used the fact that 


_9. 

9i|j 
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Introducing f” = P^^ROTj (y)g’^ and 

= p'^]p , we finally have 

X = X^+ - f2T26^)PX’^ - 

f = + P^(T^6^+T2d^-fiT^6^)Pf’^- (12) 

T -u 

- P ROT^(y)(’> T.6.)g^ 

i=l 


The Partial Derivatives of Lookpoint Coordinates 


Henceforth, we will use a prime to denote the matrix of PD with respect to d 
computed at the nominal point. From (10) it follows that 


M 


9 u 

3d 


-1 —1 -1 1 - 
= u = X + hf + h f 


(13) 


Introducing e^^ = e 2 = 1 and e^ = ^e^e ’ ^*1* be rewritten as 


’’5 e. (X.+hf.) = a 

^ — 1 1 X 

1=1 


or, 


■( 5“ f .ef) + 2h(T“ f.X.ei^) + T*X? e? = 

/C-. 1 1 XXX *— X X 


Differentiating the last expression as a inpliclt function of h gives 


h' = 


(Xi + hfh (X^- + hf-;)e? 

r f.(X, + hf.) e? 

XX x' X 

2 .„ 1 „ .L 


= _ 3 T Ui-ei(Xj+hf-}) 


f .u.e^ 
X X 1 


1 , 


and, after substitution h in (13), we have 

r 

l/k 




fl, ^ (X 

ii^k 


} + hfbu.e^l 

X X X X J 


where 


(i,k = 1,2,3) 
Co = ~y fiUje? 
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Using the matrix notations, u may be expressed as 
= C(X^ + hf^) , 


(14) 


where C is a 3x3 matrix with the elements 
C, 


, = 1 + U, f, ef C ^ 
kk k k k o 


2 ...1 

C, . = u.f e.C„ 
kj 3 k j o 


(15) 


Transformation to LSOM Coordinates 

To transform the lookpoint coordinates_j_ u, to the LSOM coordinates, they 
must be represented in the normalized form V = u* ^ . 


Differentiation of V Yields 
-2 


“U| 


’k- ” 

the elements 

\k - '-"k 


-1 


and introducing the matrix B with 


(16) 




V may be written as 


= Hull ^ BU^ = Ijulj ^ BC(X^ + hfS 


(17) 


The next step is transformation of V to W and then to X^. From (9) and (17) 


it follows that 


= 


-^1 


ABC(X^ + hfS . 


(18) 


From (7) and (8) it follows that 

1 „ 1 + Wi 

X 1 = %R*ln -j ppL- 

ml 1 - W, 


(19) 


X - = R’arctan (-W„/W„) 

m2 2 3 


and therefore, 


1 2-11 

X T = ,R(1-WJ W, 
ml 11 


4 ■ - « 3 " 2 >- 


e-10 



Introducing the matrix 


D = / 1 0 0 

0 -Wj 


( 20 ) 


and using (18), we have 


—1 2 - 1—1 

X = R(l-Wf) W = 
m i 


R ||u II ~^(1-W^) ^ DABC(X^ + hf^) 


( 21 ) 


obtain the final result, we must substitude an explicit expression for 
X + hf , which follows immediately from (12) : 


+ ■! 


'“-hp\0T^(Y) \ g" 

k= 1,2,3 

P^Tj^_ 3P(1P + hf’^) 

T — n 

= P T, „PU k = 

k-3 

—n 

-X 

k = 

— n 


k = 


Description of the Algorithm 

Calculation of the partial derivatives is performed simultaneously with the 
LSOM coordinate generation in the following order. 


Compute matrices C,B, and D, given by (15), (16), and (20) 
Compute matrices 

J = — — -rr- s DABC 


3. 

4. 


IWll (1-wJ ) 

oT 

Compute vector Z = Pu’^ 

Form five vectors 


J = JP^ 
o 


'^l r ^s2®2 ^s3®3’ ^s2^3’*'^s3®2^ 
^2 " — ^®3^’®l^s2’®l^s3^ 


“^3 - ^®2’^’®l^s3’ ®l^s2^ 

r 


= (0,-Z^,Z2) 
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= (U2,-U^,0) 

Here r = || X . 

5. Compute 

\ (k = 1 , 2 ,..., 5) 



|i y JvJ y 

Here Pj^ is the 2x1 matrix of the partial derivatives with respect to 5 ^ and thus, 
P~ ( pj^,P 2 jP 2 >P^>P 5 > Pg > Py ) • 


Note, that SCD/GCD calculations require only the first six pair of the PD. The 
partial derivatives with respect to tlme,P 2 , will be used only to generate the 
backward scan grid for Thematic Mapper. 

The nominal SCD and PD are computed in double precision and stored in single 
precision. Because the PD are changing very slowly over a frame, they may be 
computed on a sparse grid followed by linear Interpolation onto a finer grid. 

For instance, implementation of our technique for MSS requires calculation of PD 
on a 3x5 grid. 


THE NOMINAL COORDINATES AND THE PARTIAL DERIVATIVES 
FOR BACKWARD SCANS OF TM 

It should be remembered that application of the developed technique to 

Thematic Mapper data requires two sets of the nominal SCD and PD - for forward and 

backward scans. But actually only one set must be obtained by the direct lookpoint 

calculation: LSOM coordinates for, say, forward scans may be easily converted to 

LSOM coordinates for backward scans. Our calculations show also that, for sensor's 

o 

misalignments less than .1 , the derivatives are practically same for both grids; 
for bigger misalignment the second set of the derivatives can be obtained by the 
linear interpolation of the first one. 

Let X (tj^) and X (t 2 ) be LSOM coordinates for adjacent forward and backward 
scans at time tj^ and respectively. Note, that for the TM 

6 t = t„ - t, ^ 2 T ^ + T . = .132205 sec 
2 1 act round 

So, we will neglect changes in the attitude angles during 5t. 
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Considering Ax^ — X (t„)— X (t^^) as a function of changes in the spacecraft position, 
sensor pointing, ana effects of the earth rotation, we may represent it as 


+ (-TTh, \ * 

4 i 1 1 2 1 

+ 6t - (p^ + „pst + + pj 


A A 

Here 0 and 0 are fictitious pitch and roll angles, reflecting a difference in 
sensor s pointing at t^^ and t^, and is the average orbital rate during 6t. 


Here we will denote the nominal pointing vector g’^ at moments of time t] 
q, respectively. The angle between their projections onto the 
(? 2 »? 3 ) plane, ( 0 ,P 2 ,P 2 ) and (Ojq^jq^)^ can be written as 


and tr 


A 

cos 0 = 

p 


^ 2^2 ^. 3*^3 




/ 2 ^ 2,% 
(q .2 + 03 ) 


or, choosing the appropriate sign, 

^3^2 ■ 


A 

0 


A 

sin 0 = 

P 


Pjij 


(l-Pp *(l-qp' 


Analogously, 0^ may be expressed as the angle between projections of P and q onto 
the plane: 


A ^ A 
0 - sin 0 = 


^ 1^3 ^ 3*^1 

O 9 1 

(l-Pp^(l-q^)' 


For zero sensor's misalignments 

. 2 P 2 P 3 




(1-P?)^ 


A 

0=0 

t 

CONVERSION TO BASIC MAP PROJECTIONS 


It should be remembered, that completely corrected imagery eventually must be 
presented in two basic map projections, SOM and either UTM or PS, To generate 
correction data in a basic map projection, it is required to invert LSOM coordinates 
to geodetic latitude and longitude and then perform the standard mapping into 
desirable projection. 
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Noting, that the normalized look point vector, V, can be expressed through the 
geocentric latitude, (|), and the longitude, \ as 

= cosX cos(j) 

= sinX cos(|) 

= sinc^ , 


and, employing well known formula 


A. 

for the geodetic latitude ^ 


2 -2 

tan jf) = a^ tan ^ , 

one can obtain 


X 

<l> 


= arctan (V 2 Vj^^) 

= arctan | a^ b ^ V„(l-V^) 

e 3 3 J 




For a given X^, V is computed by the inverted formulae (10) and (9). 


NUMERICAL RESULTS 
The Accuracy of the Method 

To evaluate the methods accuracy, differences between LSOM coordinates, computed 
directly on points of a grid, and those, corrected according to Eq.(l), were 
calculated for various spacecraft position and attitude deviations. It is convenient 
to characterize the upper level of errors by the maximal along-track (AT) and 
cross-track (CT) errors, which coorespond to the errors at the worst points of a 
frame (possible different for AT and CT errors). It should be noticed, that the 
maximal errors always appear near the corner points and similar for TM and MSS. 

They are linearly dependent upon magnitude of deviations and practically independent 
upon WRS latitude. 

The actual position and attitude departures for Landsat-D are expected to be 
01°(a) for the pitch, roll, and yaw and less than 5km in the along and cross track 
directions. The radial departure is determined chiefly by the orbit fluctuations 
and it will not exceed 9.5km. Modeling shows, that for ^ i “ d 2 “ d ^ “ *03 , 

5 , = = 5km, and 5, = 9.5km, the corresponding maximal CT and AT errors have the 

order of 5m (CT = 4.87m, AT = 5.03m for MSS and CT = 4.97m, AT = 5.07m for TM) . 

The inversion from LSOM to geodetic coordinates produces insignificant additional 
errors, therby preserving the same order of errors in UTM and PS projections. 

For TM, the forward to backward scan conversion results in CT and AT errors . 
less than .03m for zero sensor's misalignments; for extremely large misalignments 
of the order of .5 , the maximal CT errors increase up to .5m. 
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Currently the SCD/GCD generation accuracy for Landsat-D are defined in terms 
of the^ average mean-squared errors (Im for TM and 1.5m for MSS). Table 1 represents 
the 90/^ maximal errors and the standard deviations of the average errors for Thematic 
Mapper, obtained by stochastic modeling. Here the attitude angles were normally 
distributed with zero means and o = .01 . Two cases of radial deviations were 
considered: a constant equal to 9.5 km, and a more plausible value from a uniform 

distribution (-9.5, 9.5) km. Because the errors do not depend significantly upon 
distribution of the cross and along track deviations, the latter were kept constant 
at 5 km. 400 samples were used to establish results for each case. The table 
also represents a case when PD were computed on a 3x7 grid and then Interpolated to 
a finer grid. The nominal SCD and PD for backward scans were recomputed from the 
data for forward scans. Note, that in all cases the standard deviations of the 
average errors are less 1 m and thus, the nominal SCD and PD, precomputed for mean 
orbit radius at the frame center, provide the geometric correction with the required 
accuracy. 


Timing 


On the VAX, the direct lookpoint calculations take about 11 msec per grid 
point, interpolation of PD - 1 msec, the nominal SCD to SCD/GCD correction- less 
than .5 msec, and inversion from LSOM to geodetic coordinates - 1.1 msec per point. 
Application of our technique for MSS requires interpolation PD to a finer grid, 
two corrections in LSOM coordinates, and Inversion to geodetic coordinates; 
altogether it takes about 3.1 msec per point. The direct on-line SCD and GCD 
calculation takes about 22 msec per point. 

It should be noted, that mapping to the SOM requires about 15 msec per point, 
which is considered excessive for on-line processing. This time may be significantly 
reduced if we take into consideration the fact that the LSOM closely approximates 
true SOM distances between points within each frame. The errors of the approximation 
are relatively small (less than 5m) and sufficiently regular to permit linear 
interpolation LSOM to SOM coordinates. It may be done by using a 9x9 grid of 
corrections, precomputed and stored in the Data Base (Ref. 1). 


CONCLUSIONS 

The SCD/GCD calculation technique is based on presentation of image 
distortions as a sum of nominal distortions and linear effects, caused by variation 
of the spacecraft position and attitude variables from their nominals. The 
implementation requires generation and storage of the nominal SCD and twelve (for 
MSS), or fourteen (for TM) matrices of PD for each distinct latitude of WRS, along 
one path. The maximal errors of the method do not exceed 5.1m at the worst point 
of a frame. The standard deviations of the average errors are less than Im. 

The speed of the processing and the accuracy that is achieved by this technique 
makes it an elegant solution in the production environment. 
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Table 1. 


The 90% maximal errors and the standard deviations of the average errors for 
constant and uniformally distributed radial deviations. 


Distribution 
of the radial 
deviation 

Interpolation 
of PD 

Forward sea 

ms 

Backward scans 

90% max 
errors (m) 

STD 

(m) 

90% max 
errors (m) 


CT 

AT 

CT 

AT 

CT 

AT 

CT 

AT 

constant 

no 

2.76 

2.37 

.61 

.49 

2.78 

2.36 

.61 

.49 

yes 

3.09 

2.69 

.63 

.67 

3.11 

2.70 

.70 

.72 

uniform 

no 

1.31 

1.78 

.29 

.32 

1.33 

1.78 

.58 

.58 

yes 

2.19 

1.41 

.42 

.35 

2.19 

1.40 

.53 

.48 
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The MSS Control Point Location Error Filter 
for Landsat-D.* 


I. Levine 

General Electric Company Space Division 
4701 Forbes Blvd. , Lanham, MD 20801 


ABSTRACT 

The theory and results of modeling for the MSS Control Point Location Error 
Filter are presented. The filter produces the maximum-likelihood estimates for 
average values of the spacecraft position and attitude errors during a single 
scene. The quality of the filter performance is characterized by the maximal 
cross and along-track residual errors for which probability distributions can 
be calculated analytically for a given pattern of control points. The filter 
with an automatic selection of the best set of estimates provides geodetic cor- 
rection at 90% of pixels with residual errors less than 40m for four or more 
control points and the mean-squared measurement errors of the order of 20-25m. 
The same accuracy can be preserved for eight or more control points and measure- 
ment errors of 30-35m. 


INTRODUCTION 

The ground control points (CP) , whose locations are measured on systematically 
corrected imagery and whose true coordinates are known from maps, give highly 
precise information on image displacements at each of the CP’s. The differences 
between true and measured locations provide the input to a filter, which produces 
refined estimates of the spacecraft ephemerls and attitude errors. Then these 
estimates are used for geodetic correction. 

The MSS filter theory, represented in Section I, is based on 

1) presentation of image distortions, expressed in Local Space Oblique Mercator 
coordinates, as a linear function of deviations in spacecraft attitudes and 
position (Ref.l), and 

2) recognition of the fact, that MSS processing is limited to a single scene with 

no more than 20 CP's. It is unlikely that any filter can assess the true time 
dependence in the deviations during a single scene. But we still believe that 
in some cases the MSS filter will be able to produce an reliable estimate of 
average rates. I, 


* Work performed under National Aeronautics and Space Administration Contract 
No. NAS5-25300. 
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An accuracy of these estimates is discussed in Section 2. The covariance 
analysis of the estimate errors shows that image distortions caused by the roll 
and cross-track deviations are so similar that their origins can be determined 
only by near perfect measurements. Thus, the filter is unable to produce an 
reliable estimate of both deviations. At the same time, the filter can provide an 
equivalent estimate for either variable, say, roll, which compensates distortions 
due to both sources. Analogously, for the pitch and along-track deviations. 

The analysis of covariances shows also, that initial uncertainties in rates may 
be reduced only for the equivalent roll + cross-track and pitch + along-track 
rates if there are more than 15 CP's and the mean-squared measurement errors are 
of the order of 10-15m. So, in cases of few CP's, that is of interest to us, 
only four estimates should be taken into consideration: for the yaw, radial and 

equivalent pitch and roll deviations. 

Section 3 introduces three global characteristics of filter performance : the 
maximal cross and along-track residual errors, together with combined error in 
distance. These characteristics can be obtained analytically and they establish 
upper levels of errors for any given configuration of CP's. The final formulae 
for probability distributions are presented? more details may be found in Ref. 2. 

It is known that pattern, which CP's form on imagery, have a strong effect on 
filter performance. Examples, given in Section 4, show that one of the most im- 
portant simple characteristics of CP's distribution is the maximal cross-track 
separation, which has been defined as the maximum of the cross-track distances 
between every pair of CP's. 

The examples demonstrate also, that for every pattern of CP's, measurement 
errors, and initial uncertainties in deviations, there is an optimal set of 
estimates, minimizing the residual errors. An approximate algorithm, providing 
the automatic selection of such a set, is described in Section 5. Results of 
modeling indicate that the MSS filter with the automatic selection provides the 
90% average errors less than .5 pixel (40m) for 4 or more CP's and the mean- 
squared measurement errors of the order of 20-25m, or for 8 or more CP's and 
measurement errors of 30-35m. 


I. THE MSS FILTER EQUATIONS 

Ref. II] shows that the local SOM coordinates of a frame point, X = (Xj^,X 2 ), 
may be . represented as 

X = X^ + p 6 , (1) 

where X^ is true coordinates of the point, 5 = ( 6 ^^, 62,..., 5^) is a vector of the 
spacecraft position and attitude deviations, and y is a_(2xl) matrix of the partial 
derivatives (PD), computed at the same point. Now let Z = (Zj^,Z 2 ) be the coordinates 
of a CP obtained from a map. We will assume that 

Z = X^ + 5 , (2) 

where 5 =(5j^, ? 2 ^ is a vector of Gaussian measurement errors with zero expec- 
tation values and the covariance matrix R. 
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From (1) and (2), it follows that the measured at a CP displacement, 

AX = X - "Z = pT - r , (3) 

is also normally distributed with 

E(^) = ecu's) - EC?) = u^ 
cov(AX) = EC"^^) = R 

Thus, the conditional probability density for AX can be written as 

P = (AX/6) = const'exp j^-%(AX + u6)'^R ^(AX + ud)J 

Let us assume that 6 is constant during a scene. Then a joint density of 
displacements at N control point is 

. N 

— -1 2 N n Tf 

P(AX ,AX ,...,AX /6) = const •// exp(-% (AX - 

k=l 

-u’^'S)^ r"^(AX^ - u^'s)) , 

where upper index k indicates AX and y associated with the k-th CP. 

— A 

It is known, that the maximum likelihood estimate of 6 (we will denote it as 6) 
is a solution of equation 

VL (T) = 0 , (4) 

where N 

L(^) = In P = -% (‘ax’^ - r'^CAx’^ - u^'S), (5) 

k=l 

and differential operator V is defined in Appendix. 

It is known also that 

E(6) = y (6) 

and in our case (Gaussian conditional density) the covariance matrix of 6 is 

T 

cov(5) = E((6-6)(^- &) ) = - (vL(6)v^) ^ (7) 

Note from here the summation index k will be omitted. Eq. (4) and (5) yield 

VL = -% ("M - y^) = 

= % r“^(^X - u^ ) = 0 (8) 
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and thus, the solution 6 can be written as 


A _i 
6 = M Y 
o o 


( 9 ) 


where 


M 


V T -1 

= R \X 


sr T 


Y = > R " AX 

o * — 


Noting that 


T V T -1 
VLV =-2.V R h=-M 


we also have from (7 ) that 


cov(6) = M 


-1 


.2 _2 


( 10 ) 


Now let ^2 I’® independent with the dispersions that case 

2 

R = ' 


and (9) and (10) yield 


A 

6 = M Y 

A 2—1 

cov(6) = M 


( 11 ) 

( 12 ) 


where M, Y are matrices with elements 

2 


m. 


ij 


^(y^i y^j + y^. P2j > 

^2 


(13) 


^i = 


5^ (Pli AX + ^ P2^ AX 2 ) 

^ '■ ft 


(14) 


-1 


,-l 


Elements of the matrix M will be denoted as m. , , i.e. M = (m.,). 

ij 13 


Once 6 is determined, it can be used for geodetic correction. With geodetically 

— 


corrected coordinates of a point being X + y6 , the residual error at the point, 
e = (Ej^, ^ 2 ^’ written as 

7=X+y'6-X =X + yt- (X + y^) = 
o 


= y(6-6) 


(15) 
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Thus, e is normally distributed with zero mean and the covariance matrix 
cov(e) = E |y(5-6) ” 

^ rp , T 2 


= y E {(6-6) (6-6)^ J- y'^ = M ^y^ 


( 16 ) 


Eq (16) defines local two-dimensional distribution of the ^residual ^errors at 
given frame point. It can be used also for detection of outliers , i.e. 
measurements at CP's. Erom (3) it follows, that after geodetic correction, the 
measured displacement at the k-th CP's, , can be expressed as 

e, = AX -y6 = y(6-6) +5 
k , 

As a sum of two independent Gaussian variables, also is Gaussian with zero 
mean and the covariance matrix. 


cov(e^) = Oj^yll h +R-Q 


(17) 


The two-dimensional probability density for is represented by the countour 
ellipses 

- _T -1 - ^ 2 

$(e) = e Q e = const = A 

It is well-known, that the probability that the 'point' is inside the countour 
ellipse is so the k-th CP should be treated as an outlier if 

$(i, ) > , 

K. 2 

where corresponds to a chosen confidence level. For instance, X = 9.21 for 
the 99% confidence level. 

All derived above formulae can be easily generalized to include the case when 
6 is a slow-changing function of time. Introducing the average deviations and rates 

during a scene, 

a =(aj^,a 2 ,... a^) and 3 = ( 3 ^, 62 »»-- 3^) > 

we can approximate the deviation at time t as 
6 = a + 3t 

k . 


(18) 


Now the displacement at the k-th CP at time t is 


AX^ = y^( a + 3t^) + ? 
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and the maximum likelihood function of a, $ can be written as 

- - — k k - - k -1 

L(a,g) = (AX - y(a + 3t))R 

k=l 


.,T;k k-— -;r^k. ,T 

(AX - y (a - 3t )) . 


From the above, one can obtain that the estimates of a and 3 (denoted as 
are given by 


A 

a 

= M. ^ 

Y 

A 

1 

e 

1 

Y 


(19) 


where Y is given by (14) and components of Y are 

y^- “l +4->‘21 ''V'= 




(20) 


The matrix M consists of four submatrices 

1 , ^ 


Ml = 


M M 


m' m' ' 


where M is defined by (13) and elements of m' and M* ’ , m'. , and m', ! , are 

, IJ IJ 

' r- / ""l 

“ij ^9.)t (21) 


f f 

m. . = 

IJ 


Ij 

g2 

Yi 1^9-, )t^ 


*2i "2j 


a- 2i "2j 


( 22 ) 


K — A — 

We have also that E(a) = a , R(3) = 3 , and 

yA A. 2 -1 

cov(a,3) = 


(23) 


A A A 

Introducing the estimate of deviation at time t, 6 = a + 3t, we have that 

A A A — — 

E(6) = E(a) + t • E(3) = a + 3t = 6 


and the covariance matrix of 6 is 



cov(^) = E((a-a)(a - a)'^) +2 t E((a-a)(3-g)'^) 
+ E((6 -B ) (6 -S )^) = 

„ t 2 ' ' 

= C + 2tC + t C , 

where C, c' and c' are the feil submatrices of M, V 

I f i i 



( 24 ) 


(25) 


Further it will be considered that the filter can estimate, at the most, the 
along-track (AT), cross-track (CT) , and radial (RAD) position deviations and 
rates, together with deviations and rates in the pitch (P) , roll (R) and yaw (Y) . 
and will correspond to the cross-track and along-track measurements. 


II. COVARIANCE ANALYSIS 

The covariance matrices, cov(§) and cov (a,^), completely characterize an 
accuracy of estimates, which can be achieved by the filter for a giveii configuration 
of CP's and the mean-squared measurement errors and a^. It is well-known, 
that a pattern, which CP's form on imagery, has a strong effect on filter performance, 
especially for a small number of CP's. At the same time, our calculations show 
that for M - 10 elements of the covariance matrices insignificantly depend upon a 
distribution of CP's. For aj^«#o„_^he standard deviation of the estimates are 
approximately proportional to At the present time, and are not expected 

to be less than 10 and 12m, respectively; the MSS filter will be processing up to 
20 CP's per scene. 

Table I shows the standard deviations of ^ computed for cr. = 10, = 12m and 

50 (randomly located in a frame) CP's. Comparison of the standard deviations with 
initial uncertainties in the spacecraft position and orientation, given in Table 2, 
demonstrates complete inefficiency of the filter in that case. The reason is simple; 
PD with respect to the R and CT deviations, as well as PD with respect to P and AT, 
are almost linearly dependent. As a result, the matrix M is nearly singular, and 
thus Eq. (9) can not give a reliable value of In other words, distortions, caused 
by the R and CT (or P and AT) deviations, are so similar that the difference would be 
revealed only in near perfect Imagery by near perfect measurements. 

It prompts not to estimate CT and AP deviations at all, considering the cor- 
responding image distoritons as a result of additional fictitious deviations in R 
and P, respectively. Thus, the filter should be treated as a source of appropriate 
geodetic corrections, rather than true estimates. 

The covariance analysis for time-dependent deviations shows that the filter 
is unable to produce reliable estimates of the Y and RAD rates even for N = 50: 

the standard deviations of the estimates are 4-5 times as much as their initial 
uncertainties. 
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At the same time, the filter provides mediocre estimates for combined R + CT 
and P + AT rates when w = 15 - 20 and , 02 are of the order of 10- 15m. Table 3 
shows such an example for N = 20, a. = lO, and 02 ~ 12m. Note, that initial 
uncertainties in the R + CT and P +'^ AT rates are'^.83 and .82 prad/sec (these values 
have been computed by data from Table 2) . 

Despite the fact, that in some favorable conditions the filter can cope with 
these rates, such a case will not be considered below. Being interested chiefly 
in the case of few CP’s we will take into account only estimates of the P,R,Y and 
RAD deviations. 


III. THE MAXIMAL RESIDUAL ERRORS 

The current requirements to geometric correction accuracy are specified in 
terms of .5 pixel 90% of the time. Accordingly, we will evaluate the filter 
performance by the 90% guantile of probability distribution, computed for the 
residual errors which were observed at points of some, say 15x15, grid for randomly 
distributed deviations, measurement errors, and possibly, CP’s locations. Three 
types of the 90% errors may be introduced on two-dimensional grids: for the CT and 

AT components Ojf the residuals errors and for the total residual displacements 
e = (e^ + Xhe last characteristic will be referred to as DIST. 

It should be noted, that actually all these characteristics can be obtained 
only by stochastic modeling. At the same time there are two additional global 
characteristics, which can be computed relatively simply: probability distributions 

of the maximal CT and AT errors. These distributions describe errors at the 
worst frame points and thus establish the upper level of possible errors for given 
CP's. 

A 

Let us introduce the error in the,j-th estimate, A. = ~ '^j' 

Eq, (15) may be rewritten as 

° “2j b (1 - 1,2, 3, 4 ) 


It is known that all PD increase towards the corner points of a frame. The CT and 
AT errors also reach maximal magnitudes at a corner point 111 , although it is 

never known beforehand at what specific point. At the corner points only four PD^ 
namely, P21’l'l2’^23’ *^14’ 1^®^® significant values. Moreover ,^with an error 

less than"^. 1%, they may be replaced with their maximum values, y (retaining, of 
course, correct sign). Thus, at the corner points 


m 


m 


^12 ^2 ^14 ^4 


^2 1*21 ^1 1^23 ^3 


Noting that y 2 ^ ^nd PD with respect to Y and RAD, have opposite signs at the 

ends of every scan line, we always can find a corner point, where V -<2 '^2^*^^1*14 *^4 
have the same sign (analogously, for h 2 i^l ^’^‘1 1*23 3^ ‘ indifferent to signs 

of and E 2 > we finally have that the maximal (absolute) CT and AT errors, 

Y2» 
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^ 1 “ 1 ^^12 '^2 I 1 ^14 ^4 1 
■^2 " I ^^21 ^ ^23 ^3 I 


( 26 ) 


Here we have omitted the superscript m. 


Fortunately, A„, A, and A^, A„ are practically independent and thus these expressions 
may be used separately to^derive corresponding distribution and moments (Ref. 2). 


The following are the final formulae for mean value and variance of Y . : 

./T 


E(y.) = 


7 


Var(y^) = (1 
4 

+ - S. S • 

IT 12 


2 2 2 

■ f) (S^ + sp + 

IT 1 2 

( p . arccos (r) + r - 1) , 


(27) 


(28) 


where 


r = (1 - p )'^ 

Sf = “22 


-1 
"22 

o2 2 -1 

"2 = ^^2 ^14 ^4 


m 


-1 

24 


\l^ 


-I 

22 


m 


44 


> for i = 1 


J 


^2 _ 2 2 „-l 

^21 11 

u2 m-1 

2 2 ^23 33 


P = 


liSil 




-1 

"33 


> for i = 2 


The probability distribution of y. can be written as 


Pr (y . A) = 


A 2 


1 


S IT 


J exp (-4^2-) 
0 


2Si 


(2 F (B (n + p • t - n • t) + F (B • (n - p • t - n • t) ) - 1 dt , 
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where 


n = 1 


B = A 


Sj^'r 


U 

F(u) = - r exp(- t^/2) dt, 

K 2ir V 


and 


S^, p are given by (29). 


For a given value of error. A, the corresponding probability can be easily 
computed by means of standard subroutines. Modeling has shown an excellent coincidence 
of theoretical and empirical values of E(y,) and VarCy,). Smirnov’s Test also 
demonstrates sufficient coincidence of theoretical and^empirlcal distributions. 

Eq. (26) may be used also for evaluation of the maximal residual error in 
distance, which we define as 


dm = (yJ d- y2)^^ 


We could not derive an exact distribution for d . But we have noticed that empirical 
distribution of d^ are similar to Gamma-distribution with the same means and variances 
Because E ( d^) and Var ( d^) can be obtained analytically, we have decided to 
approximate distribution of by Gamma-distribution, x\rhich is written here as 


Pr(A) = 


b^r(a) 


c a-i 

j " 


exp(- ^) 


du 


where 


a 


E(d^ 

m 

Var(d^ ) 
m 


Var(d^ ) 
m 

E(d^ ) 


(30) 


(31) 


Because, y ^ and y^ are practically Independent, 
E(d^) = E(yJ) + E(y^) 

Var(d^) = Var(yJ) + Var^y^) 

In Ref. 2 It is shown that 
2 2 2 ^ 

F(y^) = + S 2 + ^1^2 ( p apccos (r) + r) 


(32) 

(33) 


(34) 
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and 


2 2 2 2 2 3 

Var(yp = 2 (S^ + Sp + + 

+ ^ (S^ + S 2 )Sj^S 2 (r + p arccos(r)) 

- ( r + p arceos(r))^ 

7TZ 


( 35 ) 


2 2 

Here and p are given by (29). 2 ^'^' (30)-(35) yield a and b, which are used 

to compute a probability Pr(A) = Pr(A ) for any A by means of a standard subroutine 
for Gamma-distribution. Smirnov's Test shows sufficient coincidence of the 
approximations and empirical distributions for d ; differences between values of 
errors for corresponding probabilities are less than 5-7%. 


IV. EXAMPLES 

Table 4 presents means, the standard deviations, and the 90% errors for the 
maximal CT and AT errors, together with the 90% errors in distance (BIST). These 

data have been obtained by modeling (M) and analytically (T) for ~ 

and Initial uncertainties given in Table 2 (except example 9, where AT = 185m, 

CT = 35m, RAD = 65m, P = R = 120 prad, and Y = 35 yrad) . 

From 300 to 600 samples have been used to establish results for each case. The 
examples correspond to five selected configurations of CP's, depicted in Fig.l. 
Table 5 describes the examples and shows the mean-squared errors of estimates. 
Examples 1-3 and 7-8 illustrate the fact that for given configuration of CP's, 
measurement errors, and initial uncertainties, there is an optimal set of estimates 
which provides minimum errors. For distribution A, that set Includes P and R for 

distribution C it includes P,R, and Y. Examples 8 and 9 demonstrate also that 

such a set depends upon initial errors in deviations. 

As we already know, the partial derivatives p„p and p^^ represent 

the main effects of the position and attitude deviations on xmage distortions. 

At the same time, there is significant distinction between Pj^^* 1*21 ^^14’ '^23' 

when the former are almost constant in a frame, the second increase their 
magnitudes along every scan line. Thus, up to the second order effects, P and R 
estimates do not depend upon position of CP's in a scene. Roughly speaking, they 
depend upon average CT and AT displacements at all CP's. On the contrary, to 
detect effects of the Y and RAD deviations, we should observe differences of these 
displacements, so the bigger the CT distances between CP's, the bigger differences 
in corresponding p^^, and ^ 23 ’ higher an accuracy of the Y and RAD estimates 

Thus, a simple but important characteristic of CP's distribution is the maximal 
cross-track separation, H, whith we define as the maximum of cross-track distances 
between every pair of CP's. 

Example 1 shows that for small cross-track separations (H=29.7 km) the Y and 
RAD estimates are absolutely insufficient (457 prad and 292 m) and, as a result, 
the residual errors are large even for very modest measurement errors. On the 
other hand, even smaller number of CP's may lead to better results if they are 
'nicely' separated (example 4 for N = 2 and H = 169 km) . Comparasion of examples 
4 and 5 shows that an along-track shift of CP's does not affect significantly an 
accuracy of results if cros-track positions are preserved, in addition, example 6 
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suggests that a shift of CP's as a whole in the cross-track direction towards to 
the frame bounds increases errors. It implies that it is always desirable to 
have CP's placed symmetrically along the track. 

Analyzing results of modeling, we have noted that the 90% CT and AT errors 
can be approximated as 

= E(y^) + 1.5 yjVar (y^) (36) 

where E(y.) and Var(Y.) are given by (27) and (28). We have no explanation of 
that fact^ but it was^verified on a large number of cases which have shown that 
an error of such an approximation usually does not exceed 5%. 


V. AUTOMATIC SELECTION OF ESTIMATES 

As it has been shown, for every pattern of CP's, measurement errors and 
initial uncertainties, there is an optimal set of estimates which reduces the 
residual errors. Consequently, the filter's performance can be Improved if it 
will automatically select an appropriate set of estimates. Our approximate 
algorithm of selection is based on the fact that Yj^ and practically 

independent, and, bigger maximal errors almost always lead to bigger average errors. 

In our specific case, the a priori known uncertainties in P + AT and R + CT are 
always bigger than errors of corresponding estimates (at least, for mean-squared 
measurement errors less than 40m). Thus, these estimates always ought to be 
included in an optimal set. Now, all we need is to compute Yj^ twice, with and 
without the use of the RAD estimate. In the second case, the standard deviation 
of the RAD estimate must be replaced with the initial mean-squared error. Analogously 
should be computed twice to determine when the Y estimates ought to be employed. 

Table 6 presents the 90% CT, AT, and DIST errors, computed on a 15x15 grid 
as a function of cr for various number of CP's.* Results for each case 

have been established by 300-500 randomly generated sets of CP, measurement errors 
and deviations. CP's were generated so that the distances between every pair of 
them were not less than 75 km for N"^ 4, 50 km for N = 5,6 and 25 km for N 7 . 
Additional restriction forbad generation of CP's on the frame borders. The 90% 

DIST errors also depicted in Fig. 2. 

As one can see, the filter provides geodetic correction with the 90% errors 
less than 40m if cr 20m and N — 4. For = 30 m only 8 or more CP's can 
guarantee that accuracy. Note these results do not Include errors due to neglected 
uncertainties in rates. These additional errors, accumulated during 15 sec (i.e. 
with respect to the frame center) can be evaluated as 8.6 m (a) in either direction. 
Being relatively small, they do not affect significantly the total errors. 

It should be pointed out, that the automatic selection only slightly reduces 
the average residual errors. At the same time, it essentially moderates errors 
in relatively rare cases of extremely bad distributions of CP's. 


* Note, that the CT and AT errors are Gaussian and thus may be described also by 
the corresponding standard deviations. 
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APPENDIX 

a - £ X 1 matrix (vector) 

-T 

a = (aj^ ,a^, . . . ,a^) - transposed vector a 
E(a) - mathematical expectation of a 

cov(a) = E((a - E(a))(a - E(a))^) - covariance matrix of a 

Var(b) - variance of b. 


V = (■ 


85 ’ ¥5~ ’ ''‘85“^ symbolic differential 


operator defined for U = (u^^ , 0 ^ , . . . ,u^) as 


T 

VU = 


9u 


If A is a £ X £ matrix. 


!j 

95, 


3U. 


96 , 


9u 
£ 

96 £ 


V(6 ^A6)V^ = 2A 
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Table 4 


The 90% maximal CT, AT, and DIST errors (in meters) 


error 


1 



. ^ 

Example 1 

f 

1 

I 

— . } 

2 

i 

^ i 

4 

^ 

i 

5 ; 6 

T 

7 ; 

8 

i 

9 1 



M 

1 

38.3 i 

9.4 i 

9.7 

12.6: 

13.3 ,57,1 

11.0^ 

5,5 

i 

10.9 1 


mean 


j 







j 



T 

38.1 

8.6 ; 

8.6 : 

12.2 

12.2 55.6 

11.0 


11.0 : 



-M 

26.0 

4.6 

4.8 

6.6' 

6.7 :40.3 

6,1 

3,2 

i 

5,7 ; 

^1 

STD 






i 



] 



T 

26.8 

4.6 

4.6 

6 . 6 

6.6 41.6 

6.1 

3.2 

6.1 i 

(CT) 

— 

r- 




....... „ . 







M 

72.9 

16.0 

15.7 

21,5 

22.4 ,113.8 

19.4 

10.0 

18.5 1 


90% 




- 





1 



T 

75.9 

14.8 

14,8 

21.2 

21.2 ,115.0 

19.3 

9.6 

19,3 i 




M 

45.8 

42.7 

30.0 

12.9 

12.2 ; 64.5 

11.8 

11.2 

7.3 ; 


mean 








1 

... 


•< 



T 

42,1 

42.0 

30.7 

12.3 

12.2 : 61.0 

10.9 

10.9 

6.5 j 



M 

30.3 

30.7 

20.4 

6,7 

6.5 44.4 

6,0 

6.2 

3,5 i 

^2 

STD 


, - 







i 



T 

29.9 

29.9 

2.0 

6 . 6 

6.6 45.7 

6,1 

6,1 

3.6 ! 

(at) 

• ■' 

• 









■1 

< 



M 

88.0 

85.9 

59.0 

22.0 

21.0 124.5 

20.0 

19.6 

12.0 


90% 












T 

84.3 

84.2 

. 

57.9 

21.2 

21.1 1126.7 

i 19.3 

j _j 

18.9 

11.3 



M 

110.0 

86.4; 

60.2 : 

27.3 

27.8 154.8 

' 24.6 

20.6 

20.3 

^m 

90% 









r- 

(DIST) 


T 

105.5 

85. 5i 

61.2 i 

26.5 

28.9 |156.7 

i 26.2 

21.5 

' 22.1 








- .,i 


. .. 
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Table 5 


Description of Examples in Table 4. 




distribu- 


mean-squared 

I 

errors 

in 

estimates 

Example 

N 

tion 

(See ELgi) 

H 

(km) 

P 

Cyrad) 

R 

(yrad) 

Y 

( yrad) 


RAD 

Cm) 

1 

3 

A 

29.7 

14 

13 

457 


292 





Table 6 


The 90% CT, AT, and DIST errors (in meters) as a function of N and 

(^2 " ""l )• 


N 

0 ^ = 10m 

= 20m 

= 30m 

= 40m 

CT 

AT 

DIST 

CT 

1 

1 AT 

DIST 

CT 

: AT 

DIST 

CT 

AT 

DIST 

1 

19.0 

51.5 

53. 1 

33.4 

59.9 

65.7 

49.5 

1 73.3 

81.6 

68. 1 

87.4 

104. 

2 

13.4 


■ 

24.2 

47.0 

49.7 

1 

34. 9‘ 57.4 

63.0 

49.0 

65.1 

74.4 

3 

' 

10.8 

20.6 

22. 1 

19.8 

38.8 

40.8 

28.4 

47.5 

52.0 

38.2 

58.5 

65.3 

4 

10.0 

16.5 

18.3 

17.6 

31.1 

33.7 

25.0 

45.9 

49.4 

33.9 

51.2 

57.3; 

5 

9.5 

14.4 

16.2 

16.5 

29.9 

32.4 

23.4 

r ' 

43.4 

46.4 

32.5 

49.8 

55.1; 

6 



■ 


26.6 

28.6: 20.8 

; 

39.6 

42.5 

28.3 

49.2 

i 

53.3 

8 

H 


14,4 

13.6 

23.5 

25.3 

7 

18. 5i 

35.1 

37.4 

24.7 

43.0 

47.1 

10 

8.3 

11.0 

13.0 

12.2 

22.5 

24.1 

17.8; 

30.6 

33.3 

i 

22. 1 

'1 

42.3 

45.2 

15 

7.8 

8.8 

9.9 

9.9 

16.3 

m 

14.2 

24.4 

26.5 

16.9 

33.4 

35.2 
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"ONBOARD UTILIZATION OF GROUND CONTROL POINTS FOR IMAGE CORRECTION" 


J. Lowrie (Martin Marietta Corporation) 


ABSTRACT 


Future remote sensing missions require real-time knowledge of the sensor 
boresight in earth fixed coordinates for calculation of image distortion co- 
efficients and control of a pointing mount for acquisition of off-nadir data. 

An analysis of inertial navigation systems reveals an inability of these 
systems to adequately solve for the sensor boresight position due to dynamic 
misalignments between the sensor coordinate frame and the gyro coordinate frame. 

A conceptual navigation system consisting of a GPS receiver, two NASA standard 
star trackers, a NASA standard gyro package, and a landmark tracker is presented. 
The landmark tracking algorithms have been developed and analyzed, and results 
show that the position of the landmark can be determined to within two tenths 
of the sensor resolution. The navigation system has been simulated, and a 
thorough error analysis has been performed. Results indicate that this combin- 
ation of sensors can continuously solve for sensor boresight position in earth 
fixed coordinates to within 15 meters. 
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"Effective Covariance Deweighting for Precision Estimation" 

by C. E. Velez, V. L. Tate, 

Computer Technology Associates 


Abstract 


The Air Force's Sunnyvale Satellite Test Center has had a 
continuing need for near real-time high precision orbit estimates 
derived from S-Band tracking in the presence of severe atmospheric 
and geopotential modeling errors. Techniques based on sequential 
estimation using dynamically derived time-correlated process 
noise models have been developed and successfully shown to improve 
state and state covariance predictability for these cases. This 
paper will present the overall approach to sequential estimation 
currently planned for the upcoming data system upgrade to the 
current Sunnyvale system. In addition, test-bed results 
utilizing actual data taken for a medium altitude (e.g. 300 nmi) 
orbiter will be shown indicating the nature and magnitude of the 
improved performance resulting from the proposed estimator. 
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Spin-Axis Attitude Estimation and Magnetometer 
Bias Determination for the AMPTE Mission* 


R.H, Thompson^ 

Naval Electronic Systems Command 
Washington, D.C. 20360 

G.F. Neal§ 

Computer Sciences Corporation 
Silver Spring, Maryland 20910 


★★ 

M.D. Shuster 

Business and Technological Systems, Inc. 
Seabrook, Maryland 20706 


Abstract 

Algorithms are developed for the determination of magnetometer 
biases and spin-axis attitude for the AMPTE mission. Numerical examples 
of the performance of the algorithm are given. 



Presented at the Flight Mechanics/Estimation Theory Symposium, NASA 
Goddard Space Flight Center, Greenbelt, Maryland, October 27-28, 1981. 

^ Physicist, Electronic Special Warfare and Space Division 

^ Technical Staff, System Sciences Division 

Staff Scientist, Research and Development Division 
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I - Introduction 


This paper describes methods for determining spin-axis attitude 
(i.e., the direction in space of the spacecraft spin axis) and 
magnetometer biases which are being investigated for ground support of 
the Active Magnetospheric Particle Tracer Explorer (AMPTE) mission. 

The AMPTE mission will consist of two spacecraft.^ The first is the 
Ion Release Module (IRM), provided by the Federal Republic of Germany, 
which will be placed in a highly elliptical orbit with apogee at 
approximately 19 Earth radii in order to release lithium tracer ions 
outside the magnetosphere. This spacecraft will be spin stabilized at a 
rate of 30 rpm. The second spacecraft is the Charge Composition Explorer 
(CCE), which will detect the tracer ions inside the magnetosphere at 
altitudes of from 300 km to 7.5 Earth radii. The CCE will be spin 
stabil ized at 10 rpm. 

Estimation of spin-axis attitude for both AMPTE spacecraft will be 
based on the measurements of the geomagnetic field and the projection of 
the Sun line on the spacecraft spin-axis, which we take nominally to be 
the symmetry axis of the spacecraft bus. 

For the purpose of this study, the attitude sensors are assumed to 
consist of a three-axis magnetometer and a Sun sensor which measures the' 
angle between the Sun line and For simplicity it is assumed 

likewise that one axis of the magnetometer is along The other 

two axes of the magnetometer define and 

The measured quantities are taken to be 

M = magnetic field vector in body coordinates 

cos 3 = S*Y., where S is the unit vector directed from 

the spacecraft to the Sun (3 is the "Sun angle"). 
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Attitude determination activities fall into two areas: 

• Determination of spin-axis attitude 

• Determination of the magnetometer biases 

Because the orbit-apogee distance for these two spacecraft is so 
great, accurate geomagnetic field data for attitude estimation is 
available only for the segment of the orbit near perigee. This is due to 
the poor accuracy of the magnetic-field model at such high altitudes 
resulting from both the small magnitude of the geomagnetic field as well 
as from fluctuations in the field caused by extraterrestrial phenomena. 
However, because of the large spacecraft angular momenta, it can be 
assumed for both spacecraft that the spin-axis attitude at apogee will 
not differ markedly from that at perigee of the same orbit. 

Algorithms for spin-axis attitude and magnetometer bias 
determination are now being investigated. These are: 

• attitude-independent estimation of three-axis 
magnetometer biases and 

estimation of spin-axis attitude from measurements 
of the Sun and geomagnetic field angle. 

Each of these algorithms are batch estimators utilizing a long segment of 
magnetometer and Sun data. The algorithms are developed in succeeding 
sections and then tested using simulated AMPTE data. 

II - Magnetometer Bias Determination 

The attitude of the spacecraft is usually not known before the 
magnetometer biases must be determined. Here an algorithm is developed 
which determines the magnetometer bias vector by minimizing a loss 
function which is independent of the attitude. 
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The quantities used throughout this section are defined as follows: 

H.(i) = jth component of the model magnetic field in the 
*1 

geocentric inertial (GCI) system at time i 

M.(i) = jth magnetometer reading at time i 

jth component of the magnetometer bias vector, which 
is taken to be independent of the spacecraft 
position 

For the ith point, an error i5(i) is defined by the following equation: 

^(i) = 


B. 

J 


The objective of this equation is to minimize the quantity <5(i) by 
adjusting the bias vector^ to its optimal value. Thus, the loss 
function to be minimized is given by 


L(B) 


.1^ “(i)p(i) 


(2) 


where w(i) is the weight associated with the ith data point. The weights 
are assumed to be normalized to unity, that is. 


N 

I 


i=l 


0){i) 


= 1 


( 3 ) 


Determining the minimum value of L(B) first requires that its 
derivatives with respect to the components of the bias vector be set 
equal to zero: 

= 0 m=l ,2 ,3 (4 ) 

m 
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where 


^ = -4 I o)(i) 
i=l 


m 


- |B - H(i)| 


(b„ - (^) 

m m' ' 


Combining Eqs. (3-5) leads to the following results: 


I G . B, = b + F (B) 
mk k m m'--' 


(6a) 


or in matrix form. 


G B = b +\F(B) 


(6b) 


where 


'^k “ *mk«|Ji|^> - - 2 <W 


(7a) 


' <(|Hp - |i!|^)V 




m 


(7b) 

(7c) 


The bracket denotes the weighted average 


N 

<A> = I <^(i)A(i) 
i=l 


(8) 


6 , is the Kronecker delta defined as unity when m=k and zero 
mk 

otherwise. 


Eq. (6) can be solved directly to obtain the best value for the bias 
vector B. 
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General Description of the Iterative Solution 

Eq. (6) is nonlinear in jB and must be solved iteratively. The 
zero-th order (trial) solution to Eq. (6), is obtained by dropping the 
nonlinear terms in comparison to the linear terms. This approximation is 
valid only when the bias is small in comparison with the actual magnetic 
field. This point is not critical, as the iteration scheme constructs an 
accurate solution even when the trial solution is not close to the true 
solution. This will be discussed in more detail in the treatment of the 
numerical example. 

The trial solution is given by 

= G'^b (9) 

where G"^ = inverse of the matrix G 

B^^^ = trial solution 

This solution may be iterated as 

The iteration continues until 



where e = some arbitrarily small value depending on the accuracy 
desired. 
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Numerical Examples 


The AMPTE engineering data simulator^ was used to generate biased 
magnetometer data for the purpose of investigating the convergence 
properties of the iterative solution. Two cases were considered: 

B/H « 1 


and 


B/H » 1 

The first case considered was B/H « 1; in this case, 200 data points 

were used in the calculation. Data at the perigee point, at which the 

magnetic field attains its maximum value, was included. The magnetic 

field can be resolved into a component along the AMPTE spin axis, H||, and 

a component perpendicular to the spin axis, Hj^. The maximum or perigee 

values for these components are H^^^ = 240 milligauss (mO) and 
MAX 

H|| = 90 mG. The input biases were chosen to be 5 mG, 10 mG, and 15 mG 

along the x, y, and z axes, respectively. The results of the bias 
determination calculation are shown in Table 1 taken from Reference 3. 
Rapid convergence and very high accuracy is obtained. The trial solution 
B^^^ (iteration 0) initially was not accurate in the y component and 
needed to be iterated to obtain satisfactory results. Investigation of 
the case in which B » H used a subset of the data used in the first 
test. Here, 100 data points well outside the perigee region were used. 
For this test, H^^^ = 5 mG and H^^^ = 2 mG. As before, the input biases 
are 5 mG, 10 mG, and 15 mG. These results^ are presented in Table 2. In 
this case, convergence is very slow and incomplete. Improved convergence 
cannot necessarily be obtained by using standard Newton-Raphson 
techniques. 
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ITERATION 

NUMBER 

LOSS 

FUNCTION 

(mG) 

By (mG) 

(mG) 

0 

54621.0 

5.00288 

12.0278 

15.0213 

i 

5153.0 

4.98344 

9.38109 

14.9473 

2 

370.0 

5.00481 

10.1647 

15.0152 

3 

29.0 

4.99870 

9.95352 

14.9959 

4 

2.0 

5.00037 

10.0128 

15.0012 

5 

0.2 

4.99990 

9.99635 

14.9997 

6 

0.01 

5.00003 

10.0009 

15.0001 


Table 1 

Bias Determination Calculation for B/H < 1 


ITERATION 

NUMBER 

LOSS 

FUNCTION 

B (mG) 

X 

By (mG) 

B^ (mG) 

0 

24100.0 

1.8 

23 

5.3 

10 

1460.0 

3.7 

5.5 

11.0 

20 

501.0 

4.1 

6.1 

12.4 

30 

240.0 

4.4 

6.3 

13.1 

40 

133.0 

4.5 

63 

13.6 

50 

31 .0 

4.6 

6.6 

13.9 


Table 2 

Bias Determination Calculation for B/H > 1 
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Ill - Spin -Axis Attitude Determination 


Once the magnetometer biases have been chosen properly, data from 
the Sun sensor and the magnetometers may be used to determine the 
spin-axis attitude. It is assumed that the spin axis is not varying over 
the data interval examined. 


A 

The spin axis is denoted by The data are 


3(i) 

= measured Sun angle at time i 

i=l,... ,N^ 

M(T) 

= measured magnetic field at time i. 

i = l,... , 

l(i) 

= (true) Sun vector in GCI at time i, 
measured from the spacecraft to the sun 


H(i) 

= (true) geomagnetic field at time i. 

N|y| 

Note that 

there will be no requirement of simultaneous 

Sun-sensor and 


magnetometer data. 

The spin-axis (attitude) vector, is subject to the following 
constraint: 

i'i»i (12) 

The spin-axis vector is chosen to minimize the following loss 
function: 


N3 


(13) 


1 M 
^ 2 


- cos n(i)|^ ..1 X a-£ 
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where 


A = Lagrange multiplier chosen to satisfy the constraint 
equation 

oj^(i) = weight assigned to the jth magnetic field 
measurement 

“|v]{j) = weight assigned to the jth magnetic field 
measurement 

The quantity n is the angle between the geomagnetic field and the 
spacecraft spin axis given by 

n = cos"^(My/|M|) (14) 

The weights are normalized to unity 

Ns % 

i=l ^ i=l 


The spin-axis vector a is chosen to minimize the loss function 



The derivative of the loss function is given by 


(16) 


Ns 

= I - cos g(i)) S^(i) 

m i = l 

(17) 


+ (i‘M(i) - cos n(i)) M^(i) - Xa^. 


13-10 



The solution to Eq. (16) may now be written as: 


where 


mk 


= <S S, >e + <M M, 
m k S m k M 


b = <cos 3 S >e + <cos n M >„ 
m m S m M 


( 18 ) 


(19a) 

(19b) 


and the brackets denote weighted averages over the magnetometer and Sun 
data. That is, 


S 

<Cj>S “s(i) Cj(i) (20) 

Eq. (18) may be written in matrix notation as 

(A - XI) a = b (21) 

where I is the unit matrix. 

Attitude Solution 

A general solution to Eqs. (18) and (19) is constructed in this 
section. The solution to these equations leads to the spin axis attitude 
in the Geocentric Inertial (GCI) coordinate system. Again an iterative 
procedure is developed to construct a numerical solution to the 
equations. An approximate solution to the problem is to take X = 0, 
i.e., to relax the constraint that a be normalized to unity. Given this 
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approximation, Eq. (18) may be solved to obtain 





( 22 ) 


Note that this vector is not normalized. In practice this solution will 
be very close to having unit norm since even with X = 0, a^ is overdeter- 
mihed in general by Eq. (18). Thus, normal izing will lead to a 
very good approximation tor ^ (see Ref. 4). An exact numerical solution 
is generated by solving for X iteratively starting with a trial solution 
X = 0 and given by Eq. (22). 

Define the function f(X) by 


f(X) = j(X)*a(X) - 1 (23) 

Given the numerical value of j(X), the Newton-Raphson method is used to 
determine X. Differentiating Eq. (23) gives 

l^=2a-^ (24a) 

and 


^ = (A - XI) ^ 


(24b) 


The Newton-Raphson scheme gives 


j(j) . j(j-i) _ 




(25a) 


= (A - X^J^I)"^ h 


(25b) 
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Numerical Example 


The spacecraft orbit in this example is of the AMPTE type, and the 
Sun and magnetometer data used covered the perigee point. The data is 
perfect (uncorrupted by random error) as generated by the AMPTE 
simulator. The "true" value of the right ascension, a, and declination, 
5 , were chosen to be 

a = 159.67 deg (26a) 

<S = 0.0 deg (26b) 

The zero-order result as given by Eq. (22) was 

a = 159.55 deg (27a) 

6 = 0.073 deg (27b) 

in very good agreement. After ten iterations, the values changed only 
slightly, as expected, namely 

ot = 159.76 deg (28a) 

6 = 0.062 deg (28b) 

IV - Conclusions 

Efficient and reliable algorithms have been developed for spin-axis 
attitude and magnetometer bias determination for the AMPTE spacecraft. 
Using simulated numerical data it was demonstrated that the methods work 
well for AMPTE mission parameters. The present work does not address 
problems associated with noise, data rate, sensor misalignments and etc. 
These problems were investigated in references (3) and (5). 
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A MATHEMATICAL MODEL OF LANDSAT-D ATTITUDE DYNAMICS 

WITH INTERNAL MOTION 

S. D, Oh, G. W. Abshire, and J. M. Buckley 
Computer Sciences Corporation, Silver Spring, MD 

ABSTRACT 

An algorithm to model the effects of internal motion by the 
solar array and the high-gain antenna on the attitude of the 
Landsat-D spacecraft is presented here. The relative torgue 
and angular momenta arising from the internal motions are 
assumed to be attitude- independent but are considered to be 
a source of attitude perturbations. The equation of motion 
for the three-body problem is derived and then compared with 
the one-body case. The effect of the internal motion on the 
control of the spacecraft is shown in a computer study of 
the problem. 

1. INTRODUCTION 

The paper presents algorithms for modeling the effects of 
internally moving parts on the attitude of the Landsat-D 
(LSD) spacecraft. The internal motions considered here in- 
clude the rotations of the solar array to follow the Sun and 
the gimballed high-gain antenna to communicate with the 
Tracking and Data Relay Satellite (TDRS) (Reference 1) . The 
LSD system is treated as a rigid three-body system for de- 
scribing the equation of motion. Modeling the disturbance 
torques produced by moving appendages is very important for 
missions such as Landsat-D, which require accurate knowledge 
of the attitude and precise control of the spacecraft. 

The relative torques and angular momenta arising from the 
internal motions are considered as attitude-independent 
variables and as a source of attitude perturbations. The 
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external disturbance torques and the angular momenta caused 
by the internal motions are generated in a profile program 
(called PROFILE) on an IBM S/360-95 computer, where null 
attitudes are assumed and are transmitted to a truth model 
on a DEC PDP-11/70 computer that simulates the effects on 
the attitude. 

In this discussion, nonstandard rotations such as a 45-de- 
gree slew of the solar array to avoid interference with the 
antenna and the switching motion of the antenna from one 
TDRS to another are neglected. In addition to the rota- 
tional motions of the solar array and the antenna, the LSD 
spacecraft contains moving parts such as the thematic mapper 
and multispectral scanner (Reference 2). However, these 
motions are disregarded here because the motions are oscil- 
latory with a high frequency (-1 Hertz) and because they 
generate zero average angular momenta. 

Section 2 discusses the mathematical derivations of the 
equation of motion and pertinent terms such as the moment of 
inertia (MOI) tensor and the center of mass (CM) . When pos- 
sible, these terms are compared with the form for the one- 
body system used by the Multimission Modular Spacecraft 
(MMS)/Solar Maximum Mission (SMM) spacecraft. Section 3 
provides simulation results to compare the three-body and 
one-body cases. Conclusions resulting from the study are 
presented in Section 4. 

2. ANALYTICAL CONSIDERATIONS 

This section presents the mathematical modeling to describe 
the dynamic effects of the moving parts on the motion of the 
spacecraft. The equation of motion for the LSD mission is 
referenced at the CM of the entire system but is represented 
in a coordinate system that is fixed in the main vehicle. 

The CM of the entire system is calculated as a function of 
time. The MOI tensors for the moving parts are reevaluated 
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with respect to a set of time-independent axes parallel to a 
set in the main vehicle. Also calculated are the angular 
velocity of the appendages and the perturbation in the ex- 
ternal torques due to the changing positions of the append- 
ages. A comparison with the one-body problem is made. 

2.1 COORDINATE SYSTEMS AND TRANSFORMATIONS MATRICES 

The system under consideration, shown in Figure 1, consists 
of the main carrier vehicle, designated as body B^, and 
n(=2) moving bodies Bj (j=l, n) . Several coordinate sys- 
tems are convenient for discussing the relative motions. 
These are as follows: 

• Geocentric Inertial Coordinate System (GCI) (Refer- 
ence 3) 

• Orbit-Defined Coordinate System (OCS) where X 
(roll) is nearly along the spacecraft velocity vec- 
tor, Y (pitch) is along the orbit normal vector, 
and Z (yaw) is along the nadir vector 

• Spacecraft-Fixed Coordinate System (BCS) , which is 
fixed in the main vehicle B 

o 

• Coordinate systems fixed in moving parts such as in 
the solar array (SACS) or in the high-gain antenna 
(ANTCS) 
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NOTE : Qj = the CM of 

T = the CM of the entire system 
CM 

■q^ = the CM of Bj from r^^^^ 

~ the hinge point of B^ 

iT. = the angular velocity of B^ in inertial space 

■£“. = the angular velocity of B . relative to the main 
^ body Bq (oTj = Wg +1?^) 

Figure 1. Partitioning of the Satellite Into Main Body 
and Moving Parts 
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The transformation matrices (TRMA) to be used in this paper 
are defined as follows; 


1. TRMA from GCI to PCS : [0] 


(Rj X 


I> ®I 



X 



I 


II 


'Cs 




K 

X 



I 


I 




^ 1 


R-r 

X 

Vx 


I 


1 1 






R-r 




I 

— 


( 2 - 1 ) 


where and denote the spacecraft position relative to 
the Earth and velocity unit vectors in the GCI frame, re- 
spectively . 


2 . Attitude direction cosine matrix from the PCS to 
the BCS : [A] . In the PROFILE Program [A] is given by the 

identity matrix because null attitudes are assumed. In the 
truth model, it is represented as 



( 2 - 2 ) 


using the small angle approximation, which is sufficient and 
valid, since only small perturbations are assumed; r, p, and 
y denote roll, pitch, and yaw angles in radian units, re- 
spectively . 
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fSAl 

3. TRMA from BCS to SACS : [C ] . The solar array 

rotates around the y-axis and is driven to follow the Sun. 
Thus, its orientation is determined from the Sunline angle a 


[C 


(SA) j 



cos a 0 . -sin a 

0 10 
sin a 0 cos a 


(2-3) 


^ T . , 

Given the Sun unit vector, S = (S , S , S ) , in the BCS, 

X Y Z 

the rotation angle a is given by 



(2-4) 


because the Sun vector is perpendicular to the x-axis of the 
SACS. 


I ant) 

4. TRMA from BCS to ANTCS ; . [C ] . The antenna 

has two gimbals with the inner gimbal angle, g^, repre- 
senting the elevation angle and the outer gimbal angle, 
g^ , representing the azimuth angle. The orientation of 
the antenna is determined from the gimbal angles 


[c(ANT)] ^ 


^11 

-I y L J Z 


COS g^ cos g^ 
-sin g^ 


sin g^ cos 
cos g^ 


-sin g 2 
0 


:2-5) 


cos g^ sin g 2 


sin g^ sin g 2 


cos g 2 J 


The unit vector pointing from the spacecraft to TORS is re- 
presented by P where P = (P^, P,,, P„)'^ in the BCS. 

X Y z 
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The gimbal angles are thus given by 


g^ = tan ^ (P /p ) 
1 y X 


and 



. -1 
-sin 


(2-6a) 


(2-6b) 


Since g^, g^ should align the antenna boresight (the x-axis 
in ANTCS) with the normalized pointing vector (P can be 

obtained from the spacecraft and TDRS ephemer ides . ) 

2 . 2 ANGULAR VELOCITY OF MOVING PARTS 

The angular velocity of the moving parts is used to calcu- 
late the internal angular momentum of the spacecraft for use 
in the equation of motion. It is easily seen from Equa- 
tion (2-3) that the angular velocity of the solar array is 
as follows; 




(2-7a) 


The time derivative of the rotation angle a can be com- 
puted numerically 


da _ g(t) - g(t - At) 
dt At 


(2-7b) 
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Using Equation (2-5) the angular velocity of the high-gain 
antenna is 



(2-8a) 


where 


(t) - g^ (t - At) 
dt " ~S€ 


(2-8b) 


For SMM, the angular velocity of the moving parts was not 
calculated . 


2.3 CENTER OF MASS 

For LSD, the CM of appendage in the BCS is given by 


Qj (t) 




X . 
3 


+ X. 
3 


(2-9) 


where Q.„ represents the CM of B. at the initial time (see 
3O 3 

Figure 1). The rotation (or hinge) point is denoted by X^ 

^ ^ J 

and QjQ - Xj represents the CM of B^ from the hinge point 
at the initial time. Then, at any later time, the CM will 
be represented by the first term of the right-hand side of 
Equation (2-9). The CM of each appendage changes as a 
function of time because the high-gain antenna rotates to 
track the TDRS , and the solar array rotates to track the 
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changes in 


Sun. Consequently, the CM of the system, 
time and is represented by 


n 





(2-lOa 

r = 0 / 

r = 0 



and the position of the CM of each appendage with respect to 
the CM of the system is 


qj(t) = Q.(t) - r^j^(t) 


(2~10b) 


For SMM, the CM of the system was fixed in time in the BCS. 

2.4 MOMENT OF INERTIA TENSOR OF THE SYSTEM 

The MOI of the system, [I^], relative to axes parallel to 
the BCS axes passing through r^^^ is expressed by 


r=0 ' 

- [q,(t) + [q^(t) -f- p^]^j 


( 2 - 11 ) 


where is the position vector of the mass dm of body 
B relative to the CM of B and the subscripts 1 and m 
represent the 1 and m components of the vector or tensor. 
Note that because is time-dependent, is also de- 

pendent on time; in the remainder of this paper, the ex- 
plicit time-dependence will be dropped. 
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The above equation can be written as. 
n 

ln\ ~ l^r 

r=0 ' 

since 



[I ] is the MOI tensor of represented in the BCS frame 
but relative to the CM of B^: 


q^(5 , - (q ) 1 (q ) 

^r Im '^r^ 1 ^r m 


+ 

y( r ) 

_ 

L- J 


J Imj 


( 2 - 12 ) 


[I = [C ^ [I [C (2-13) 

/ y- \ 

where [I^ '] is the MOI of B^ represented in the coordinate 
system fixed in B^. . Equation ’(2-12) can be simply reexpressed 
by 




(2-14) 


r = 0 


with 


(r)' 


Im 


= + M !q^5. - (q ) , (q ) i (2-15) 

L Jim r P^r Im ^r 1 ^r m ) 


For the one-body problem, as represented by SMM, I is de- 
fined to be a constant in time. 

2.5 EXTERNAL TORQUES 

TWO external torques are discussed: the gravity gradient 

torque and the aerodynamic torque. The solar radiation 
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torque is similar to the aerodynamic torque, and the other 
external torques are not sensitive to the three-body problem. 

The gravity gradient torque, can be computed by 





R + p + q 

r ^r , 

X — : — - dm 

I R + + "q i 

I ' r . ' 


(2-16) 


where y is the Earth gravitational constant (==3. 986005 x 

10^“^ m^/sec^). R is the spacecraft position vector from the 

Earth. Considering that l^l>> 1^^ + ^l, is simply, 

r r CjyC:t 



The expression for the one-body system has the ‘same form 
except for the replacement of [I^] by the constant [I]. 

To simplify the calculation of the solar radiation and aero- 
dynamic torques, the LSD spacecraft is modeled as an as- 
sembly of a cylinder for the main vehicle, flat plates for 
the solar array panels, and a sphere for the antenna. Only 
the aerodynamic torque is discussed here because the modi- 
fications to the center of pressure (CP) are common in solar 
radiation and aerodynamic torques. 
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The aerodynamic torque, is 



( 2 - 18 ) 


Here, v denotes the spacecraft velocity unit vector, n^ 
denotes the normal unit vector for the ith surface, • 

O ^ J. 

denotes the CP of the ith surface from p denotes 

the atmospheric density, and Cj^ denotes the drag coeffi- 
cient. The normal vectors, n^, for the solar array and 
antenna surfaces are dependent on time by 


n . 



T 




lO 


( 2 - 19 ) 


where n. 

10 

surface . 
puted by 


represents the initial normal vector for the ith 

'q for the solar array and antenna are com- 

^cp,i 


with 


^cp,i ®cp,i ^CM 


( 2 - 20 ) 


‘cp 




(i) 


cp , 10 


Xi) 


+ X. 


( 2 - 21 ) 


More consideration is required to specify for the 

solar array surfaces that are canted. The transformation 
matrix from BCS to these surfaces, is given by 


[C 


= [C 


(SA) 




with the canted angle 9 ^. 


( 2 - 22 ) 
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For the one-body case of SMM, 
2.6 EQUATION OF MOTION 


and ^ are constants. 


The equation of motion for the LSD spacecraft is written in 
the form 


dt 


f (Y(t) , t) 


(2-23) 


where Y = (q^, (U = If 2, 3, 4) denotes the 

Euler symmetric parameters representing a rotation from the 

GCI to the spacecraft-fixed coordinate frame, iT^ is the 

total angular momentum of the spacecraft, and “l,, is the 

w 

wheel momentum. 

The body angular momentum of the main vehicle, is 

B 

given by the total spacecraft angular momentum minus the sum 
of the wheel momentum, payload momentum, ^ , and the angu- 
lar momentum, caused by the internal motions 



(2-24) 


Lg depends on the angular velocity of the main vehicle, 
o)q, and depends on the angular velocity of moving 

parts. To"!. To formulate these mathematically, the 
angular momentum of the total system, L^, ignoring wheel 
and payload momenta, is considered 


^T S ^INT 


n 


/ ^^r^ ^ ("qj. X ^^) dm^ 


r=0 

n 


= Ur X q'j. + (oj^ + oj;.) 


o r ' 


r=0 


(2-25) 
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with some computation, can be shown as 


^ = 


[1^] % ^ 

r=0 


where 


[K^ l^r * ^'^CM “ ^r^ ^Im " ^ ’^CM " ^r^^ ‘^m 


Thus, the body rate of the main carrier is simply 


a>o = 


and caused by the internal motion, is 

INT 


'INT 


. r = l 




Ul 


The time derivatives of the Euler symmetric parameter 
q , can be obtained as 

y 


dq 


dt "2 ^ yv 


with 


[ m ] 


1 

0 

“3 

-‘^2 

I 

rM 

3 

-0)3 

0 

<"1 

0)2 

“2 

1 

£ 

!-■ 

0 

^3 

_-0)i 

CM 

3 

1 

■ -‘^3 : 

£ 

1 


( 2 - 26 ) 


( 2 - 27 ) 


( 2 - 28 ) 


( 2 - 29 ) 


( 2 - 30 ) 
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The time derivative of the total angular momentum of the 
spacecraft is given by the Euler equation as 


dt ~ ^ext ^ ^o 


( 2 - 31 ) 


For SMM, the body angular momenta, L-, is given by 

D 




( 2 - 32 ) 


with the payload momentum L„. The spacecraft body rate, 

K 

1^, is determined by 


/Ui-, 


0) = I (i)2 I = [I' 

,0J. 


B 


( 2 - 33 : 


where [I]~^ is the inverse of the spacecraft MOI tensor. 
The time derivatives of the Euler symmetric parameters, 
q , can be obtained as 


dq 


dt 2 ' ' HV 


: 2 - 34 ) 


The time derivatives of the total angular momentum of the 
spacecraft are given by the Euler equation as 


— = “ext + I-T " 


( 2 - 35 ) 


with the external torque, N 


exf 
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SIMULATION RESULTS 


A computer study of the effect of the three-body problem on 
the motion of the spacecraft has been made using the general 
equations derived here. Since the spacecraft is subject to 
noticeable external torques, a control law that provides 
compensatory torques was necessary to keep the spacecraft 
near null attitude. The one-body case, using the same con- 
trol law, was also studied. 

The roll, pitch, and yaw of the spacecraft main carrier for 
both cases is shown in Figures 2 through 4. The results of 
the three-body case are represented by the "X" points and 
the results of the one-body case are shown as open circles. 
Note that both cases are subject to the same control law. 
This control law attempts to make the pitch, roll, and yaw 
zero and to bring the spacecraft rate to null. This control 
law is the same one (Reference 4) that Landsat-D will use 
during its acquisition phases. The torque applied to each 
reaction wheel is as follows: 

for the roll axis. 


T 


r 


K (k Ar + 0 ) ) 
r r r 


for the pitch axis, 

Tp = Kp [kp(4p + B) + 


and for the yaw axis. 




(3-la) 


(3-lb) 


(3-lc) 


where Ar and Ap are the roll and pitch attitude errors as 

determined by an Earth sensor: K,K,K,k,k,k, and k 
i 'r p y r p y 
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Figure 2. Spacecraft Roll Versus Time 
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are constants; B is a bias to compensate for the orbital ro- 
tation; and (jj^, o)p, and are the angular velocity along 
the roll, pitch, and yaw data. Because of the values used for 

k , k , and k , the control law is much more sensitive to 
r P Y 

the spacecraft rate than to the attitude error. 

Most of the structure seen in the plots is a result of the 
control law. However, since the control law is the same, 
the differences in the plots are a result of the three-body 
problem. Note in Figure 2 that after 4.5 minutes the con- 
trol law has the roll rate to zero for the one-body problem 
but not the three-body problem. Likewise, after 2.5 min- 
utes, the pitch rate of the one-body problem is under con- 
trol. 

4. CONCLUSIONS 

The conversion of the rigid one-body problem to the three- 
body problem has added another dimension to the study of 
dynamics. Although the exact perturbations in motion are 
obscured by the control law used, the effects are still im- 
portant in control of the spacecraft. 

The algorithms used in this paper can be applied to other 
spacecraft such as the Space Telescope to study important 
low-frequency effects, as in this paper, and also higher 
frequency effects that will cause jitter. 
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SPACECRAFT ATTITUDE POINTING PERFORMANCE DURING ORBIT 
ADJUST AS A FUNCTION OF COMPENSATOR ORDER 


For many communication satellite missions, it is required that the 
control system performance during velocity adjust mode does not degrade 
appreciably from the nominal pointing requirement. During velocity adjust, 
many factors contribute to the development of disturbance torques that exceed 
the capacity of the reaction wheels. This necessitates the use of thrusters 
to provide the control torques. The spacecraft weight constraints force the 
use of off-pulsing techniques. While off-pulsing the orbit adjust thrusters 
may eliminate propellant penalties, it also introduces additional disturbances. 
The thruster plume impingement torques increase dramatically when the bal- 
ancing effect of both thrusters firing is lost. 

In order to meet the attitude pointing error requirements under a set of 
constraints outlined above, a steady state compensator of specified order is 
proposed to estimate the required duty cycle needed to balance the disturbance 
torque. The compensator order has been increased gradually to demonstrate the 
improvement in pointing accuracy. The basic mathematical model of the flexi- 
ble spacecraft and sensor used to characterize the performance of the compen- 
sator can be described as follows: 


e(t) = (02 + 02 )h - 2^0302 ?(t) - 0)^02 <P(t) 

(1) 

* 

H(t) = T - T 
d c 

(2) 

0 (t)= - 1 0 (t) + 1 0(t) 
Ti ^ 

(3) 

• 

i 0,(t) + G 0. (t) 

(4) 

<p(t) = c(.t) 

(5) 

h(t) = -25coi;(t) - ui^ip(t) + 02H(t) 

(6) 

Td = 0 

(7) 

y(t) = 02(t) + V (t) 

(8) 
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where , 



5 

CO 

? 

y 

V 

G 


sensor time constant 

disturbance torque 

control torque 

spacecraft attitude 

sensor output after first break 

sensor output after second break 

the spacecraft momentum 

- 1/2 

spacecraft rigid body admittance = (Inertia) 

structural admittance at first symmetric (pitch) or 

asymmetric (roll/yaw) frequency 

structural damping 

structural frequency 

modal deflection 

integral of modal deflection 

noise corrupted sensor measurement 

measurement noise 

sensor gain 


The continuous model of the estimator has been represented as 


a = Fg, + h y 

(9) 

T 

u = -g z 

(10) 


where matrices F, h define the compensator structure and g is the feedback 
gain. The vectors ^ and ii define the compensator state and the control 
respectively. 
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The problem presented in this paper involves estimating the distur- 
bance torque using a compensator of specified order as represented in 
equations (9) - (10). As a baseline, the compensator is assumed to be a 
third order to estimate the rigid body position, the momentum and the dis- 
turbance torque. The compensator order is gradually increased to estimate 
the sensor states and the flexible modes. Having specified the dimension 
of the compensator, the matrices F, h and g have been chosen to minimize 
the performance criterion involving quadratic function 



The performance criterion for this problem has been chosen as 

J = Lim E(L) 
t-1^ 

where E (* ) denotes expection. 

The attitude pointing performance has been documented as a function of 
the dimension of the compensator. The analysis thus provides a trade-off 
between increased pointing accuracy and increased complexity in on-board 
software. 
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Efficient Algorithms for Single -Ax is 
Attitude Estimation^ 


M.D. Shuster^ 

Business and Technological Systems, Inc, 
10210 Greenbelt Road 
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Abstract 

Computationally efficient algorithms are presented for determining 
single-axis attitude from the measurement of arc lengths and dihedral 
angles. The dependence of these algorithms on the solution of trigono- 
metric equations has been much reduced. Both single-time and batch 
estimators are presented along with the covariance analysis of each 
algorithm, 

^Presented at the Flight Mechanics/Estimation Theory Symposium, NASA 
Goddard Space Flight Center, Greenbelt, Maryland, October 27-28, 1981. 

^Staff Scientist, Research and Development Division 
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I. Introduction 


Since nearly every spacecraft is spinning during part of its 
life — in particular, at the time of orbit injection — spin-axis attitude* 
estimation is an important segment of almost every mission support 
operation. Indeed, for spin-stabilized spacecraft there is often no need 
(or desire) to determine the complete three-axis attitude at every 
point and, in fact, when accuracy requirements for the spin-axis attitude 
dictate that many measurements taken at different times be processed 
simultaneously, the computation of a three-axis attitude may not even be 
possible. 

Very often, three-axis attitude information is definitive data 
required chiefly by mission scientists and generally processed anytime 
from several days to several months after the receipt of telemetry. The 
need for efficient three-axis attitude estimation algorithms in those 
cases is determined by the definitive data rate. When three-axis 
attitude information is required in real-time for the purpose of attitude 
control, this is usually provided on-board by three-axis gyros (e.g. SMM) 
or on the ground by the spin axis and a third angle, which can be 
obtained by monitoring some other sensor reading such as IR scanner pitch 
(e.g. AEM, Magsat). 

Spin-axis attitudes by contrast are usually required not only as 
definitive data but also by the ground support system in near real-time 
for the purpose of monitoring spacecraft performance and determining 
large scale attitude maneuvers. Thus, the efficiency of a spin-axis 
attitude estimation algorithm becomes a factor in the safety and daily 
operation of the spacecraft. 


* 

Since the single-axis attitude of interest is invariably the spin-axis 
attitude these terms will be used almost interchangeably throughout this 
work. 
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While a number of highly-efficient algorithms exist for three-axis 
attitude estimation,^ the computation of spin-axis attitude^ is by 
comparison very clumsy. This is largely because the computation of 
three-axis attitude uses complete vector measurements in general and can 
take advantage of the linear properties of Euclidean three-space. The 
computation of spin-axis attitude, on the other hand, must rely on 
incomplete vector information (the measurement of arc lengths and 
dihedral angles) to determine a quantity (the spin-axis) which is 
restricted to the surface of a sphere. Thus, while three-axis attitude 
computations need only execute simple matrix operations, the computation 
of spin-axis attitude is beset with the burden of solving complex 
relations from spherical trigonometry. 

Since spin-axis attitude is usually not computed frequently, the 
need for efficient algorithms is not immediate, at least not for ground 
support systems. The determination of the spin -axis attitude from batch 
measurements of arc lengths and dihedral angles has become highly 
standardized and reliable^ and there is no obvious need to replace this 
software in normal ground support operations. 

The need for more efficient algorithms lies in two areas: 1) the 

eventual implementation of spin -axis attitude computation in onboard 
microprocessor-based attitude determination systems; and 2) the computa- 
tion of spin -axis attitude accuracies, which imposes a far greater 
computational burden than computing just the attitude due to the greater 
number of terms and because the computation of the attitude covariance 
involves implicitly the computation of derivatives of the attitude. 

The large computational burden imposed by the need to solve 
spherical trigonometric equations in the computation of spin-axis 
attitude covariances is evident in the work of Wertz and Chen,^**^"® 
the most complete and careful work to date. The difficulties which are 
encountered in this approach are of two kinds: 1) the complexity of the 

trigonometric relations, themselves, and 2) the fact that for certain 
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cases the representation of the quantities being calculated becomes 
indeterminant while the quantities themselves are well defined. This 
last difficulty is simply a manifestation of the fact that the 
representation of rotations by Euler angles is sometimes ambiguous and is 
overcome in the same way, namely, by changing the representation. 

The need for computing spin-axis attitude covariance matrices is 
two-fold. Firstly, it is necessary to be able to assess the accuracy of 
a spin-axis attitude computation during the spacecraft mission. 

Secondly, it is important to be able to predict spin-axis attitude 
accuracies for mission planning, particularly in the determination of 
launch windows. For an example of launch window computations using the 
geometrical approach see Chen.^ 

The purpose of the present work is to develop algorithms for 
computing spin-axis attitude and the associated covariance matrix without 
relying as heavily as do current methods on the solution of trigonometric 
equations. A completely vectorial approach is, of course, not possible 
owing to the nature of the measurements themselves. However, in large 
degree many of the trigonometric equations can be abandoned with the 
result that the spin-axis attitude and, particularly, the covariance 
matrix can be computed more efficiently. 

The types of measurements studied here are of two kinds: 

measurements of arc length , which will always be the angle 
between the observed direction and the spin axis. 

measurements of dihedral angles , i.e., the angle between two 
planes, where the line of intersection is assumed to be the 
spin axis.® 

Dihedral angles, in general, are measured by observing two crossing 
times in the spacecraft and multiplying by the angular velocity. Arc 
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lengths may be measured in a variety of ways, for example, by direct 
sighting (as Of the Sun or a star) or by measuring the component of a 
vector along the spin axis (e.g., the magnetic field vector). The 
measurement of the nadir angle is hybrid in that an arc length (the nadir 
angle) is determined from the measurement of a dihedral angle (the Earth 
width). It is the measurement of the nadir angle which is the source of 
most of the computational complexity. 

Estimation algorithms may be classified either as deterministic 
(usually single-frame, i.e., single-time) algorithms, in which a minimal 
subset of the available data is chosen to compute the spin-axis attitude, 
or as optimal (batch) algorithms, in which a larger quantity of data is 
used from which one computes a "best" result. Three cases are treated in 
this report 

1) A deterministic estimator using two arc -length measurements, 

2) A deterministic estimator using the measurements of two arc 
lengths and the included dihedral angle. (Since in this 
case the spin-axis attitude is over-determined the question 
of optimality is also discussed.) 

3) An optimal batch estimator utilizing any number of 
measurements of dihedral angles and arc lengths. 

In each case the covariance analysis is presented in detail. 

In the appendix the measurement of the nadir angle is presented. It 
is at this point that trigonometric relations cannot be avoided, at least 
in so far as measuring instruments (horizon scanners) are presently 
constructed. The treatment is similar to that of Wertz and his 
collaborators (Ref. 2) but a method is given for avoiding sign 
ambiguities. 
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The treatment of single-axis attitude estimation presented here 
complements that of Wertz. The advantage of Wertz's treatment is that 
the variances along two great circles of the celestrial sphere 
intersecting at the direction of the spin axis and the dihedral angle 
between these two circles (the correlation angle) is given fairly 
directly. Much less direct is determining the covariance of the 
spin-axis vector in inertial space. This part of the calculation falls 
out simply in the present formalism. 

The results presented here are quite simple although they do not 
seem to be generally known. An important result, which is demonstrated 
here, is that little accuracy is lost by relaxing the constraint in the 
optimization that the spin -axis vector be a unit vector and then 
unitizing post hoc. This is responsible for a great deal of 
simplification of the methods presented here, especially for batch 
estimation. 


II. Single-Frame Spin-Axis Estimation from 
the Measurement of Two Arc Lengths 

Consider the simplest case in which the measured quantities are e, 
the Sun angle (the angle between the spin axis and the Sun vector), and 
n, the nadir angle (the direction between the spin axis and the nadir 
vector). The case where one of these measurements is replaced by the 
magnetic field angle is analogous. 

Let ^denote the Sun unit vector, E the nadir vector, and ^ the spin 
axis. Then 

S*n = cos 0 = c- (la) 

E*n = cos n = c_ (lb) 

W bM# I- ' ' 


16-6 



The direction of the spin-axis can then be determined simply by using a 
method that has been published recently by Grubin,^ though it has been in 
use since the beginning of the space program and probably has been known 
for several hundred years. 

If ^and J^are not parallel, then it is always possible to write 



i = aJ + a^E + a„ i X E 

(2) 

The problem is 

now to determine the coefficients a^, a^, a^^. 


From Eqs, 

(1) and the normalization condition we have 



C3 = H = + a^d-D 

( 3 a) 


= id= 85 ^E 

( 3 b) 


1 = i.-a = 3s + + 2a3a3(S*E) + a^|i x J|2 

( 3 c) 


which have the solution 




li^i 


2 


■ I *. ^ ?.|2 


S X E 




S X E 


- (Cs - 2C5C^(H) + 


(4a) 


(4b) 


(4c) 
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Note that there are two possible solutions for n. These are shown 
geometrically in Figure 1. 

It will be convenient to define the following quantities 





( 6 ) 


where the tilde below the letter denotes a two-dimensional ‘vector or a 
2x2 matrix, 

Eqs, (4) can now be written 

a = U c ( 7a ) 



S X E 


[1 


- c'^U 


(7b) 


The covariance analysis is now straightforward. Define the three 
vector 



( 8 ) 



Then the covariance matrix of the measurements is given by 

s <6c 6c^> (9) 

where the bracket denotes the expectation value and 6c is the error in c. 
The covariance matrix of the spin-axis direction in the non-orthogonal 
coordinate system is 

Pg = <«a (10) 

and in an orthogonal coordinate system 

( 11 ) 

Substitution of Eqs. (7) in Eq. (10) gives readily 



wi th 

M = <6a 6a^> = U P^ 

~ -c ~ 

V = M b 
S = b^ M b 
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(13a) 

(13b) 

(13c) 
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where 


B = 


I 


I I b 

I ~ 



(18) 


Equations (17) and (14) may now be combined to give 


where 



2 

I 


i=l j=l 




Xi = 

w-1 


I.+ 


(19) 


(20a) 



» E + 


I) 


Eq. (16) 1s again satisfied since 


(20b) 


,X.-n = 0 i=l,2 (21) 

III. Single-Frame Spin-Axis Estimation from the Measurement 
of Two Arc Lengths and the Included Dihedral Angle 

The ambiguity in determining the spin-axis observed in the previous 
section is removed if the included dihedral angle is also measured. The 
dihedral angle t is defined as the angle between the (S,n) and (£,n) 
planes and is easily shown to be given by 


sin i|> 




■'(i-(5-i)^)(i-(£-a)^; 


(22a) 
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(22b) 


cos = 


( 1 * 1 ) - 

Al-{^n)'^)(l-(E•n)^) 


tan ij> = 


A A A 

Ji* (1 £) 


(H) - {H){|*i) 


(22c) 


The geometry is depicted in Figure 2. 


To determine the spin axis attitude it will be convenient to define 


=-i{l-Cs){l-c|) sin ♦ 


(23) 


and 


c = 


(24) 


The vector a is now determined by four equations 




(25a) 




(25b) 


A 

= s 


X eI^ a 


(25c) 


1 = * a| * 2 ajajCS-l) + ag 


(25d) 
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The three components of ^ are now overdetermined. The most convenient 
solution is obtained by solving the first three equations, which are 
linear, leading to 

a = U c (26) 

WWW 

where 


U = 


S X E 


/N A 

1 0 

-(H) 1 0 

0 0 1 


(27) 


The spin-axis in given by this however, is not properly normalized 
since the measurements are not exact, A properly normalized spin-axis 
vector is then obtained by simply normalizing the solution 


"=ii/|a| 

The covariance matrix of a is given simply by 


(28) 


P = U 
a c 


(29) 


and the covariance matrix for the unnormal ized spin -axis is given by 


P„.TP^tT (30) 


similarly to Eq. (14). The covariance matrix of the properly normalized 
spin-axis vector is recovered simply as 



QPnQ 


(31) 
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where 


Q = I - n (32) 

It is well to ask how good is the approximation of ignoring the 
normalization condition and then normalizing the solution post hoc. 
Instead of this seemingly brutal approach one can find the best solution 
to Eqs. (25abc) subject to the constraint of Eq. (25d), i.e., one seeks 
to minimize the loss function 


L(a) = (c-Aa)"^ (c-Aa) 


subject to the constaint 


where 


a^ A a = 


1 



1 (S*E) 0 

(S-i) 1 0 


(33) 


(34) 


(35) 


The solution is straightforward and yields 

a„pt - (A c (36) 

where X is the Lagrange multiplier for the constraint and from Eq. (34) 
is the root of the equation 

c c 

which yields the smallest value of the loss function. 
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Equation (36) may be rewritten 


Jopt ° (I-lPcU)-'a 


(38) 


where a is given by Eq. (26). Since is expected to be close to a, 
it follows that XP^U must be small. An approximate solution for a^p^ can 
be obtained by expanding Eqs. (37) and (38) in xp and solving. This 
yields 


Now 


-opt 


1 p n a 

2 aVA-^a “ 

c 


(39) 


<l-a^Aa> = Tr(Pj.U) 


(40a) 


<(l-a^Aa)^> = 4 aVa 

' MM' «nf ^ mM 


(40b) 


so that the additional root mean square (rms) error in ^ when optimality 
is not taken into account is of the same order of magnitude as the rms 
error in the cosine measurements. However, the source of this additional 
error, as shown by Eqs. (40) is the error in the normalization. Hence 
this error will be greatly reduced when the unit vector is normalized. 

IV. Batch Estimation 

The value of avoiding trigonometric expressions becomes more obvious 
in dealing with batch estimation. The computational advantage of the 
present approach over the geometrical approach^ is substantial. 
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For batch estimation the non -orthogonal basis cannot be used since 
only the Sun vector is constant (and then only for relatively short data 
spans). The present treatment focuses on the case where the measurements 
consist of two arc lengths and the included dihedral angle. The 
extension to other cases is straightforward. 


Let Cj(i), Cg(i), Cj^(i) be a series of measurements of the Sun 
projection, the nadir projection, and the Sun-nadir dihedral angle, 
respectively. Then the best solution for the spin-axis is obtained by 
minimizing 



subject to the constraint 



( 42 ) 


In order to decrease the number of subscripts in the expressions it has 
been assumed that each data type is available at each time and that each 
measurement type has a single characteristic error. Except for a 
proliferation of subscripts the expressions which follow are not changed 
when this assumption is removed. 
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The minimization of Eq. (41) subject to the constraint is 
straightforward and leads to 

n,= (M-Xir’v (43) 


where 


M = 


'"1 "s "e “n 


(44a) 


V = 


1=1 <>s ‘^E 


(44b) 


and X is the root of 


^1= 1 

(M-XI)2 


which leads to the smallest value of Eq. (41), 


(45) 


As in the previous section it can be expected that the constraint 
can be ignored (X»«0) and the solution be approximated by 


i = i!/|l»| (46) 

where 

n = 1 (47) 
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This approximation has been tested for one spacecraft^® and been observed 
to be quite good. The covariance of j]^ is given by 

P„ = «■' (48) 

and the covariance of the normalized solution is given again by 

P = Q P„ Q (49) 

n 


V. Measurement Errors 

The computation of the spin-axis covariance matrix requires as 
input a model for the covariance matrix of the cosine measurements. 
Expressions are derived here for computing these for the case of Sun and 
Nadir measurements. The treatment when one of these measured quantities 
is the magnetic field is treated in the same way. 

Sun Measurements 

The quantity measured is usually the Sun angle, 3. Hence, 

5Cj = -sinB6& (50) 


Nadir Measurements 


If the spacecraft has angular velocity w, then the Earth width is 
given by 


« = “(tg-tj) (51) 

where tj and t^ are the in- and out-triggering times, respectively, of 
the Earth scan (for a momentum-wheel mounted scanner, o) will be the 
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angular velocity of the momentum wheel). 


Then, using the results from the appendix 


= 6cosn 


3cosn 

H «cos^ 


9C0S^ 


Sinn ^ 

cotY - cotn ^ ? 


= ^ (sin (6t^ - 6t.) (52) 

2 cotY - cotn ' 2^ ' 0 I' ^ ' 

where y is the scan-cone half angle. 

Dihedral Angle Measurements 

The dihedral angle t is determined from the time interval from the 
Sun crossing to the mid-point of the horizon scan 

= 0)[tj - -^(tp + tj)3 (53) 

Thus, (3,Ji,'l») or (0,n,t) is a set of statistically independent 
variables. The "dihedral cosine” Cn^, however, is given by 

C|^| = sin3 Sinn sirnj^ (54) 


hence 


6C|y| = C|nj[cote 60 + cotn 6n + cott 6i|>] (55) 

From Eqs. (50-55) the covariance matrix can easily be calculated. 
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To a large degree, much of the trigonometric complexity which has 
been removed from the attitude solution has simply been shifted to the 
computation of a derived measurement covariance matrix. There is, 
however, a substantial gain because the covariance matrix need not be 
computed to the same degree of accuracy as the spin-axis attitude 
itself. Hence, a great deal of computational approximation is possible, 
such as approximation of the trigonometric functions by simple rational 
functions. 


Appendix - Measurement of the Nadir Angle 

Because the Earth is an extended body the nadir vector is not 
measured directly but determined from measurements of the Earth width. 
Earth widths are measured by a horizon scanner, which effectively is a 
sensor mounted on a rotating cone (of half-cone angle y) about the 
spacecraft spin axis, which detects the crossings of the Earth horizon on 
the scan cone. The Earth has an effective angular radius of p, which is 
a function of altitude and (for a non-spherical Earth) latitude. The 
Earth width is the dihedral angle between the in- and out-crossings (Hj 
and Hq) the horizon by the scanner and is denoted by n. These quantities 
are related by the spherical law of cosines^ 

cosp = COSY cosn + sinY sinn cos(J2/2) (A-1) 

The geometry is depicted in Figure 3. 

Eq. (A-1) may be solved to give 


cosn 


COSP COSY ± 


sinp cos 


(n/2) 4 


-cos P 


A 


(A-2) 
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where 


A = cos^p + sin^Y cos^(n/2) (A-3) 

The sign ambiguity may be eliminated if another measurement is 
present, say that of the Sun angle, 3, and the Sun-Earth dihedral angle, 
Y. Let C be the arc length from the Sun direction to the mid scan point 

cos5 = cos3 COSY + sin3 sinY cosi{> (A-4) 

Then it is possible to show that the underdetermined sign in Eq. (A-2) 
must be the same as that of 

(cos3 - cosy) (E*S - cos5) 

Alternatively, one may consider simultaneously Sun and horizon 
measurements. This leads to three simultaneous equations 

cosB cosn + sin$ sinn cost = (A-5a) 

COSY cosn + sinY sinn cos(n/2) = cosp (A-5b) 

2 2 

COS n + sin n = 1 (A-5c) 

Equation (A-2) was obtained by solving Eqs. (A-5b) and (A-5c) 
simultaneously. One could just as easily solve Eqs. (A-5a) and (A-5b) 
for cosn and sinn. The result will not necessarily satisfy Eq. (A-5c) 
but the two equations have the advantage of being linear. The solutions 
can then be renormalized to satisfy Eq. (A-5c). 

This approach of ignoring the proper normalization for the 
trigonometric functions has another advantage in that a simultaneous 
solution to Eqs. (A-5b) and (A-5c) may not exist in certain extreme cases 
because the measurements are not exact. By solving Eqs. (A-5a) and 
(A-5b) a solution will always exist. 
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There is, however, one clear disadvantage. If Eq. (A-2) is used then 
e, n, and will be statistically independent. If, however, the linear 
equations are solved, n will be correlated with 0 and Thus, the 
simplicity gained in computing cosn is counterbalanced by greater 
complexity in computing the measurement covariance matrix P . 



Figure 1 

Single-Axis Attitude from Two 
Arc-Length Measurements 
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Figure 2 

Single-Axis Attitude from Measurements 
of Two Arc Lengths and One Dihedral Angle 


A 

A 



A 

A 



Figure 3 

Geometry for Nadir-Angle Determination 
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